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ABSTRACT 


M 

In  structural  engineering  It  is  Imperative  to  design  each  system 
to  survive  the  Inputs  anticipated  over  the  design  life  of  the 
structure.  Strong  motion  inputs  cause  systems  to  execute  nonlinear 
responses,  and  during  the  strong  motion  responses,  structures 
accumulate  damage.  Therefore,  the  capability  to  model  nonlinear 
responses  and  to  assess  the  damage  level  in  a  structure  is  essential 
for  optimal  design. 

Techniques  for  the  diagnosis  of  damage  In  inelastic  structures 
have  been  developed.  The  dissipated  energy  in  mechanical  systems  is 
taken  as  a  measure  of  damage  accumulation.  Two  models  for  the 
simulation  of  damaged  structural  response  have  been  developed.  Both 
the  single-degree-of -freedom  and  mu lti-degree-of -freedom  systems  were 
Included  in  the  analysis.  The  objective  of  this  study  is  to  use  these 
models  to  estimate  the  amount  of  energy  dissipated  due  to  a  strong 
motion  input. 

The  results  show  that  structural  damage  can  be  predicted,  even  in 
the  presence  of  measurement  noise. 
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CHAPTER  1 
INTRODUCTION 


1.0  Introduction 

When  a  structure  is  excited  by  an  external  force,  it  executes  a 
response  determined  by  the  characteristics  of  both  the  input  and  the 
structure.  One  could  predict  the  exact  response  of  a  structure,  char¬ 
acterized  by  its  geometry  and  its  mechanical  properties,  if  he  could 
predict  inputs  exactly;  if  he  had  a  perfect  model  for  the  structure; 
and  if  the  mathematical  computations  were  correct.  However,  since 
inputs  are  random,  one  cannot  perfectly  characterize  complex  struc¬ 
tures,  and  since  mathematical  models  are  not  perfect,  one  can  only 
estimate  the  response  of  a  structure. 

In  structural  analysis  it  is  necessary  to  assess  the  response  of  a 
structure  to  dynamic  loads,  such  as  blasts  and  earthquakes.  This 
procedure,  of  course,  requires  the  use  of  a  dynamic  model  which  will 
permit  accurate  prediction  of  the  response  of  a  structure.  These 
structural  models  are  generally  chosen  to  fit  experimental  data -and  to 
simplify  mathemtical  computations. 

Most  existing  structures  were  designed  based  on  a  static  model, 
and  although  dynamic  properties  may  be  considered  in  their  design,  the 
design  parameters  mqy  be  Inadequate  to  predict  the  response  to  dynamic 
load  correctly.  Considerable  work  has  been  performed  on  identifying 
the  parameters  of  mathematical  models  from  dynamic  experimental  data, 
and  various  approaches  have  been  proposed  for  predicting  system  param¬ 
eters  based  on  experimental  data. 

These  identified  parameters  can  be  used  to  predict  the  dynamic 
response  of  a  structure  to  a  different  excitation  than  that  used  to 
test  It.  The  identified  parameters  also  can  be  used  to  calculate  the 
energy  dissipation  in  a  hysteretlc  structure  caused  by  strong 
excitation. 


The  energy  dissipated  by  a  structure  during  a  strong  motion 
response  is  an  indication  of  the  structural  damage.  It  is  important  to 
predict  how  much  damage  occurs  in  a  structure  due  to  strong  motion 
because  the  level  of  damage  is  related  to  the  likelihood  of  structural 
failure.  When  damage  does  occur.  It  can  appear  in  different  forms, 
such  as  cracks,  permanent  deformation,  or  change  in  characteristic 
frequency. 

In  practice,  it  may  be  difficult  to  assess  the  damage  extent  and 
location  in  a  complex  structure  after  an  extreme  excitation.  For  exam¬ 
ple,  in  a  nuclear  power  plant  or  burled  protection  structure,  it  may  be 
difficult  to  assess  damage.  This  difficulty  may  arise  due  to  the  num¬ 
ber  of  elements  in  a  structure  or  the  scale  of  individual  structure 
members. 

A  certain  degree  of  damage  is  usually  unavoidable  when  structures 
are  subjected  to  strong  motion;  therefore,  estimation  of  structural 
damage  Is  necessary  for  proper  design.  At  present,  the  exact  criteria 
useful  in  judging  the  failure  of  a  structural  system  are  not  avail¬ 
able.  It  Is  known  that  some  of  the  measures  of  structural  response 
related  to  the  occurrence  of  failure  are  peak  response,  energy  dissipa¬ 
tion,  plastic  deformation,  etc.  In  fact,  the  failure  criteria  of  any 
practical  material  are  a  complicated  function  of  many  measures  of 
response. 

In  this  study  models  are  established  for  si ngle-degree-of -freedom 
(SDF)  and  multi-degree-of-freedom  (MDF)  damageable  structures.  Signals 
that  can  be  used  to  simulate  measured  data  are  generated.  These  are 
used  to  Identify  the  parameters  of  the  structural  models.  Finally, 
damage  measures  are  computed  for  the  simulated  systems  and  their  mathe¬ 
matical  models. 

1.1  Literature  Review 

The  present  Investigation  establishes  mathematical  models  for  the 
simulation  of  damaged  structural  response.  The  mathematical  models  are 


proposed,  then  signals  which  can  be  used  to  simulate  the  input  and 
response  of  damaged  structures  are  generated.  These  signals  are  used 
to  identify  the  model  parameters.  There  are  two  broad  areas  in  the 
literature  that  are  concerned  with  concepts  important  in  this  study. 
These  are  (1)  mathematical  models  for  damaged  structural  systems,  and 
(2)  the  identification  of  structural  system  parameters.  Some  papers 
from  the  literature  in  both  these  areas  are  briefly  discussed  below. 

Part  of  the  energy  dissipated  by  a  structure  is  dissipated  due  to 
hysteretic  behavior  of  the  structural  material.  The  equation  governing 
the  hysteretic  response  of  a  lumped  mass  system  is  a  second-order, 
nonlinear,  ordinary  differential  equation  with  history-dependent  stiff¬ 
ness  term.  Two  models  which  may  approximate  the  nonlinear  system  will 
be  proposed  in  this  study.  These  are 

1.  High-order  equivalent  linear  system; 

2.  Time-varying  parameter  linear  system. 

The  first  model  considered  in  this  investigation  is  a  high-order 
equivalent  linear  system.  It  is  assumed  that  the  nonlinear  hysteretic 
system  is  approximately  governed  by  a  high-order  equivalent  system. 

This  model  is  motivated  by  studies  summarized  in  the  literature.  For 
example.  Lutes  and  Hseih  [1]  used  a  third-order  linear  system  to 
approximate  a  SDF  oscillator  with  bilinear  hysteretic  yielding  behav¬ 
ior,  excited  by  stationary  white  noise.  In  the  linear  system,  certain 
parameters  were  chosen  so  that  the  root  mean  square  displacement  and 
velocity  matched  empirical  values  for  the  nonlinear  system.  They 
showed  that  the  third-order  system  gives  a  better  overall  prediction  of 
response  buildup  than  does  either  the  linear  SDF  system  or  a  two-mode 
system. 

Lutes  [2]  used  a  different  type  of  equivalent  linear  system  to 
approximate  the  nonlinear  system.  All  the  methods  Lutes  considered 
defined  the  equivalence  either  in  terms  of  response  displacement  level, 
velocity  level,  frequency,  or  a  combination  of  these.  He  found  that 


a  particular  equivalent  linear  system  can  generally  only  be  expected 
to  match  a  limited  number  of  response  statistics  of  a  particular  non¬ 
linear  system  with  a  particular  type  of  excitation. 

Wen  [3]  and  Wen  and  Baber  [5,  6]  have  used  the  equivalent  lineari¬ 
zation  method  to  approximately  represent  the  response  of  a  hysteretic 
SDF  system.  They  showed  that  the  third  order,  linear,  differential 
equation  provided  a  satisfactory  representation  of  the  inelastic, 
hysteretic  systems.  This  closed  form  linearization  is  relatively 
simple  to  formulate  which  allows  ready  extension  to  multi-degree-of- 
freedom  (MDF)  systems.  They  showed  that  the  equivalent  linearization 
method  gives  satisfactory  results  at  all  response  levels  for  response 
analysis  of  MDF  deteriorating  or  non-deteriorating  systems  under  random 
excitation. 

Another  study  by  Wafa  [7]  demonstrated  that  the  peak  response  for 
an  hysteretic  SDF  system  excited  by  random  inputs  is  closely  predicted 
by  a  third  order,  linear  equivalent  system.  Recent  work  [8]  has  also 
shown  that  the  high-order  linear  equivalent  model  provides  a  good 
approximation  to  the  hysteretic  system  when  the  energy  dissipated  and 
frequency  shift  are  concerns.  Significantly,  the  results  established 
that  the  parameters  of  a  higher  order  system  can  be  identified  by  using 
a  frequency  domain  method  even  when  noise  is  present  both  in  the 
forcing  and  response  signals.  In  contrast,  the  time  domain  approach 
yields  poor  results  in  the  presence  of  noise.  Because  of  the  frequency 
domain's  preferable  application,  it  will  be  used  to  do  the  most 
analysis. 

The  second  model  is  motivated  by  the  fact  that  structures  may 
exhibit  time-variant  nonlinear  response  to  strong  motion.  This 
implies  that  structure  deterioration  was  in  progress  during  the  large 
amplitude  motion.  Such  a  phenomenon  has  been  recognized  and  studied  In 
the  past.  For  example,  Udwadla  and  Trlfunac  [9]  and  Iemura  and 
Jennings  [10]  carried  out  an  analysis  to  characterize  such  behavior  in 


terns  of  a  quasi-time-variant  linear  formulation  based  on  the  data 
obtained  from  the  San  Fernando  earthquake  of  February  9,  1971. 

In  another  study  Townsend  and  Hanson  [11]  demonstrate  time-varying 
hysteretic  loops  by  the  experimental  test  of  reinforced  concrete  beam- 
column  and  T-shaped  specimens  under  different  loading  conditions.  In 
addition,  Uzumeri  [12]  has  also  shown  the  same  behavior  for  an  experi¬ 
mental  study  of  cast-in-place  reinforced  concrete  beam-column  joints 
subjected  to  simulated  seismic  loading. 

Based  on  the  above  referenced  investigations  involving  time- 
varying  parameters  systems,  we  anticipate  that  the  time-varying  param¬ 
eter  model  will  provide  a  good  representation  of  a  hysteretic  system. 

Many  papers  in  the  literature  address  the  problem  of  damage 
assessment;  however,  few  of  these  treat  the  mathematical  quantification 
of  damage  measures.  In  the  following,  some  papers  which  discuss  damage 
analysis  In  both  quantitative  and  non-quantitative  ways  are  discussed. 

To  assess  structural  damage.  It  is  necessary  to  first  define  and 
quantify  the  damage.  Yao  [13]  has  examined  various  definitions  of 
structural  damage  and  reviewed  available  methods  for  damage  assess¬ 
ment.  In  1971,  Wiggins  and  Moran  [14]  developed  a  procedure  for  grad¬ 
ing  existing  buildings  In  Long  Beach,  California.  Damage  is  assessed  on 
a  point  basis,  and  a  total  of  up  to  180  points  Is  assigned  to  each 
structure  according  to  the  evaluation  of  structural  components  of  five 
types.  In  1975,  Culver  et  al  [15]  proposed  the  field  evaluation  method 
(FEM)  which  is  applicable  even  when  building  plans  are  unavailable.  In 
1980,  Bresler  et  al  [16]  described  their  structural  and  fire  evaluation 
model,  which  was  developed  to  provide  a  broad  overview  of  potential 
safety  problems  for  more  than  10,000  buildings  for  a  government  agency 
in  the  United  States.  Recently,  Ishlzuka,  Fu,  and  Yao  [17,  18,  19] 
suggested  a  rule-based  damage  assessment  system  called  SPERIL  Version 
I.  All  these  systems  are  primarily  based  on  professional  experience 
and  engineering  judgment  in  the  decision  making  process. 


A  more  mathematical  quantification  of  damage  in  structures  has 
been  used  by  Ang,  and  Wen  [20].  They  used  the  hysteretic  energy 
absorbed  and  the  maximum  structural  distortion  as  the  function  of 
structural  damage.  Others  Mho  published  in  this  area  are  Rudd,  Yang 
and  Manning  [21],  Yao  [22],  Toussi  and  Yao  [23,  24],  Chen  and  Yao  [25] 
and  Yao,  Toussi,  and  Sozen  [26]. 

An  important  aspect  of  the  present  study  is  the  method  used  to 
Identify  the  parameters  of  the  damaged  structure.  The  literature  on 
structural  identification.  In  general,  is  very  broad.  A  few  of  the 
papers  closely  related  to  the  present  investigation  are  reviewed  here. 

The  historical  development  of  research  in  the  area  of  system 
identification  is  summarized  in  the  works  of  Astrom  and  Eykhoff  [27], 
Bekey  [28],  Bowles  and  Straeter  [29]  and  Collins,  Young,  and  Kei fling 
[30].  Many  survey  papers  have  been  written.  For  example.  Col  lings,  et 
al  [30],  Sage  [31],  Rodeman  and  Yao  [32],  Chen  [33],  Hart  and  Yao  [34], 
Ting,  Chen,  and  Yao  [35]  and  Liu  and  Yao  [36]  present  surveys  of  struc¬ 
tural  identification. 

The  potential  for  change  In  structural  characteristics  due  to  the 
accumulation  of  damage  exists  and  can  be  investigated  through  observa¬ 
tion  of  structural  parameters.  Signature  analysis  technqlues  have  been 
used  to  predict  cracking  In  bridges  by  Cole  [37,  38]. 

One  Important  area  in  system  identification  permits  the  analyst  to 
characterize  system  modes.  An  early  paper  by  Kennedy  and  Pancu  [39] 
shows  how  model  parameters  can  be  obtained  from  a  vector  representation 
of  the  steady-state  response  In  the  complex  plane.  Later  papers  in  the 
same  area  were  written  by  Bert  [40],  Bert  and  Clary  [41],  Smith  [45], 
and  Trail-Nash  [46].  Some  of  the  difficulties  encountered  in  applying 
the  techniques  developed  in  the  above  papers  are  described  by  Nord 
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The  least  squares  approach  to  the  identification  of  system  param¬ 
eters  has  been  developed  in  many  studies,  and  is  applied  in  both  the 
time  and  frequency  domains.  (This  technique  is  used  in  the  present 
study.)  Some  papers  that  have  been  written  on  the  subject  of  least 
squares  parameter  estimation  are  those  of  Distefano  and  Rath  [48,  49], 
Flanelly  and  Berman  [50],  Hart  and  Yao  [34],  Ibanez  [51],  Ibrahim 
[53-58],  Milne  [59],  Raggett,  Rodeman,  and  Yao  [32],  Ting,  Chen,  and 
Yao  [35],  Udwadia  and  Shaw  [61],  and  Wells  [62].  Brieman  [63]  shows 
that  for  a  linear  time  invariant  system,  the  least  square  prediction  is 
optimal.  The  technique  for  using  the  least  squares  parameter  identifi¬ 
cation  in  the  frequency  domain  was  presented  by  Ibanez  [51,  52].  Wang, 
Paez,  and  Ju  [8,  66,  67]  found  that  frequency  domain  approach  is  well 
suited  to  the  identification  of  parameters  for  high  order  linear  models 
and  time  varying  linear  models. 

The  effects  caused  by  measurement  noise  can  be  important  in  the 
identification  of  structural  system  parameters.  Kandianis  [68]  consid¬ 
ered  that  effect  for  linear  structural  systems.  His  theoretical  devel¬ 
opment  considered  white  noise  random  input.  The  studies  by  Wang,  Paez, 
and  Ju  [8,  66,  67]  show  that  the  parameters  of  higher  order  linear  and 
time  varying  linear  models  can  be  accurately  predicted  even  in  the 
presence  of  noise. 

1.2  Objective 

The  determination  of  the  system  parameters  from  suitable  experi¬ 
mental  observations  is  a  fundamental  problem  in  engineering.  Obtaining 
a  good  representation  of  a  system  requires  all  the  proper  information, 
such  as  well  measured  data  and  a  suitable  model. 

The  objective  of  this  study  Is  to  justify  two  possible  models  to 
characterize  the  behavior  of  a  system.  The  relative  merits  for  each 
model  are  discussed.  Extensive  numerical  experimentation  using  simu¬ 
lated  data  is  also  presented  In  order  to  investigate  their  feasibility 
and  accuracy.  This  study  will  demonstrate  how  well  the  system  param¬ 
eters  can  be  identified  with  and  without  noise  in  the  measurements. 


Once  the  parameters  are  known,  the  energy  dissipated  by  the  system  can 
be  computed.  Based  on  the  computed  results,  one  can  compare  how  well 
the  models  performed  for  a  given  set  of  data.  The  ultimate  goal  of 
this  study  is  to  establish  structural  models  useful  for  other  purposes, 
such  as  prediction,  design,  control,  and  damage  assessment. 

The  present  research  has  been  aimed  at  the  analysis  of  damage 
accumulation  in  concrete  structures.  It  is  assumed  that  when  a  con¬ 
crete  structure  dissipates  energy,  it  accumulates  damage.  To  justify 
this  assumption  some  physical  experiments  have  been  performed.  Specif¬ 
ically,  concrete  cylinders  have  been  subjected  to  cyclic  loading.  The 
energy  dissipated  in  each  cylinder  was  measured  and  the  level  of  resid¬ 
ual  strength  in  each  cylinder  was  determined  after  the  load  cycling  was 
completed.  The  residual  strength  was  plotted  versus  energy  dissipated. 
When  the  reduction  in  strength  Is  taken  as  a  measure  of  damage,  this 
plot  reveals  the  damage  caused  by  energy  dissipation. 


CHAPTER  2 

HIGH  ORDER  EQUIVALENT  LINEARIZATION 


2.1  Model 

The  differential  equation  governing  the  response  of  a 
single-degree-of -freedom  (SDF)  system  is 


nff  +  u  *  f 


(2-1) 


where  m  is  the  mass  of  the  structure,  f  is  the  forcing  function, 
z  is  the  displacement  response,  dots  denote  differentiation  with 
respect  to  time,  and  u  is  the  restoring  force  of  the  structure. 
Equation  (2-1)  can  be  used  to  model  the  actual  system  in  which  u 
I  can  be  a  very  complicated  function.  In  the  present  study,  t‘  » 

hysteretic  restoring  force,  u,  is  modeled  by  using  the  equation 


W  +  2 


(2-2) 


where  the  cj,  j  *  0,  1,  ...,  M+l  are  the  constants  governing 
the  system  restoring  force  characteristics,  where  u(j)  denotes 
the  jth  time  derivative  of  u,  and  M  is  a  constant  denoting  the 
order  of  approximation  provided  by  the  linear  system.  The  reason 
for  using  this  model  to  represent  the  hysteretics  system  is  that 
it  displays  a  hysteretic  character  that  can  be  made  to  match  the 
character  of  an  inelastic  structure. 

Combine  Equations  2-1  and  2-2  in  the  following  way.  Solve 
Equation  2-1  for  u,  then  take  derivatives  of  the  resulting 
expression.  Use  these  in  Equation  2-2;  the  result  is 
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(-f{J)  +  mz^+2))  -  f  -  ui i  (2-3) 


where  fu'  is  the  jth  derivative  of  f,  and  zu  C}  is  the  (j+2)th 
derivative  of  z.  To  simplify  the  identification  procedure. 
Equation  2-3  was  divided  by  co*  This  particular  arrangement  is 
chosen  since  Co  should  never  be  0  for  the  systems  under 
consideration. 

Consider  the  case  where  M  is  equal  to  0.  The  model  in 
Equations  2.1  and  2.2  become? 

mz  +  —  2+_Lz  =  f  (2-4) 

c  o  C  o 

This  is  simply  the  second-order  linear  differential  equation 
governing  the  SDF  system.  However,  vrfien  the  response  of  the 
actual  system  is  linear  and  damping  is  viscous,  the  model  of 
Equation  (2-4)  represents  the  actual  system.  The  restoring  force 
function  for  this  system  is  u  =  (ci/co)z  +  (l/co)z.  This  model 
displays  the  hysteretic  behavior  as  shown  in  Figure  (2-1). 

When  the  constant  M  is  chosen  as  1  in  equation  2-2,  the 
model  becomes 

mz  +  u  =  f 

COU  +  ClU  =  C2Z  +  z  (2-5) 

Combining  these  equations  results  in 

rife*’  +mz+-££z+~-z*f+-£if  (2-6) 

Co  Co  Co  Co 

The  parameters  of  the  system  In  Equation  (2-6)  can  be  chosen  so 
that  the  model  represents  the  hysteretic  system  as  well  as 
possible. 

For  example.  Figure  (2-2)  shows  the  hysteretic  properties 
for  an  SOF  system  by  plotting  the  restoring  force  versus  dis¬ 
placement.  The  parameters  for  the  system  and  the  forcing  input 
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Figure  2-1.  Total  restoring  force  versus  displacement  for  second-order  system 
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are  given  in  Figures  (2-2).  This  study  will  consider  both  the 
equations  (2-4)  and  (2-6). 

The  parameters  can  be  identified  by  using  the  least-square 
identification  criterion.  Since  most  observed  data  include  a 
certain  percentage  of  noise,  the  frequency  domain  approach  will 
be  used  to  perform  the  parameter  identification. 

It  is  anticipated  that  as  the  order,  M,  of  the  model  in 
Equation  2-2  is  increased,  the  response  of  a  hysteretic  system 
can  be  matched  with  increasing  accuracy.  However,  for  practical 
reasons  involving  estimation  accuracy  for  system  parameters,  very 
high  order  linear  system  models  cannot  be  used  to  simulate 
hysteretic  system  behavior. 

.  Before  considering  the  problem  of  parameter  estimation  for 
the  system  of  Equation  2-2,  note  that  we  anticipate  calculating  a 
set  of  parameters  with  values  in  a  specific  range.  For  example, 
when  M  -  0  and  Equation  2-4  is  the  model  for  system  response,  we 
anticipate  finding  values  1/co  >  0  and  ci/co  >  0  (i.e.,  co  >  0 
and  ci  >  0.  These  values  guarantee  that  the  model  has  positive 
stiffness  and  damping,  as  we  know  the  real  system  must.  When 
M  =  1  and  Equation  2-5  models  the  system  response,  we  anticipate 
finding  values  co  >  0,  ci  >  0,  and  C2  >  0.  These  values  guaran¬ 
tee  that  the  model  response  will  be  stable. 

The  energy  dissipted  by  the  systems  described  in  this  sec¬ 
tion  can  be  computed  using  Equation  2-1  or  2-2.  In  the  terms  of 
Equation  2-1,  the  energy  dissipated  by  a  system  is 


z(T) 


where  z(0)  is  the  system  displacement  at  time  0,  z(T)  is  the 
system  displacement  at  time  T,  and  T  is  the  time  through  which 
the  energy  computation  is  executed.  This  formula  provides  the 
area  enclosed  in  the  hysteresis  loops  generated  by  the  system 


response.  For  example.  Equation  2-7  could  be  used  to  compute  the 
area  enclosed  by  the  hysteresis  loop  of  Figure  2-1.  This  equa¬ 
tion  can  be  transformed  into  terms  more  convenient  for  computa¬ 
tion.  Note  that  the  integral  of  Equation  2-7  is  written  in  terms 
of  the  displacement  variable  z.  The  variable  of  integration  can 
be  transformed  to  a  time  variable  yielding  the  following 
expression. 


(2-7a) 


Finally,  Equation  2-1  can  be  used  to  obtain 
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(f  -  mz)  z  dt 


(2— 7b) 


When  the  input  and  response  for  an  SDF  system  are  measured. 
Equation  2-7b  can  be  used  to  directly  compute  the  energy  dissi¬ 
pated  by  a  system.  When  the  response  is  computed,  for  example  by 
solving  Equations  2-5,  then  the  input  and  system  parameters  are 
used  to  find  z  and  its  derivatives.  The  input,  f,  and  computed 
response,  z,  are  used  in  Equation  2-7b  to  determine  the  energy 
dissipated. 


Numerical  examples  where  we  compute  the  parameters  of  linear 
systems  equivalent  to  bilinear,  hysteretic  systems  are  presented 
in  Chapter  8.  In  these  examples,  the  energy  dissipated  by  each 
system  is  computed  and  the  results  are  compared. 


CHAPTER  3 

IDENTIFICATION  OF  PARAMETERS  IN  THE  TIME  DOMAIN 


3.0  Time  Domain 


roach 


The  time  domain  parameter  identification  procedure  is  used 
to  introduce  the  identification  procedures  for  the  models 
described  in  Chapter  2.  Though  time  domain  approach  is  not  used 
extensively  in  this  investigation,  the  method  can  still  be  effec¬ 
tive  under  certain  circumstances.  Particularly,  the  time  domain 
identification  process  is  useful  when  the  measured  input  and 
response  contain  little  or  no  noise. 

Measured  input  and  response  data  from  a  structural  system 
are  used  to  estimate  the  system  parameters.  The  measured  input 
data  can  be  a  base  acceleration,  force  or  a  pressure  function. 

In  the  following,  the  measured  response  data  are  assumed  to  be 
given  as  acceleration  values.  This  assumption  is  realistic  since 
structure  response  acceleration  is  often  measured  during  an 
experiment  or  test  by  accelerometers  installed  in  the  structure. 

3.1  Formulation 

U) 

Consider  Equation  2-3  and  let  z^,  zA,  and  ,  £=0,...,n-l, 
be  the  response  displacement,  velocity,  and  jth  derivatives  of 
the  displacement  at  time,  t&  *  £At,  *  0,...,n-l.  Let  f£  and 

fp),  £  *  0,...,n-l,  be  the  force  at  time  t^,  £  *  0,...n-l  and 
its  jth  derivative.  Then  Equation  2-3  can  be  written  at  time 
t*  to  obtain 
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£  -  m££ 


£*0,...,n-l  (3-1) 

The  reason  for  normalizing  with  respect  to  Co  is  that  co  should 
never  be  0  for  the  system  under  consideration. 


The  notation  in  this  equation  can  be  simplified  by  taking 
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(3-2) 


Equation  3-1  can  be  written  as 

M 

a°2i  *  a‘2»  ♦  S  aj*i  I 

j-i 


(‘  fIU)  *  1 f,  '  ">*,  . 


A  «  0,...,n-l  (3-3) 


The  notation  in  this  equation  can  be  simplified  by  defining 
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(3-5) 


where  the  T  superscript  refers  to  matrix  transposition.  Using 
these  expressions  in  Equation  3-3  yields  the  relation 


M  !*l  ‘ f*  -  1 . . . . 


(3-6) 


This  Is  the  equation  governing  the  linear  system  response  at 
time  t£. 

The  notation  can  be  further  simplified  by  defining 


r  f2o>i 

tzi} 


fo  -  nUo 
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Using  this  expression,  the  sequence  of  Equation  3-6,  for 
*  -  0,...,n-l,  can  be  written 

N  M  ■  !fz!  <3-8> 

This  equation  governs  the  linear  system  response  at  all  times. 
When  (1)  the  system  from  which  the  data  were  measured  is  truly 
linear,  (2)  there  is  no  noise  in  the  measured  data,  and  (3)  the 
derivatives  and  integrals  of  f £  and  z£,  *  *  0,...,n-l,  are 
known  exactly.  Equation  3-8  can  be  satisfied  exactly  by  the 
measured  data  and  a  set  of  coefficients.  In  general,  these 
conditions  cannot  be  satisfied  exactly,  therefore,  an  error  term 
should  be  included  in  Equation  3-8.  Define  the  error  vector  as 


The  element  e£,  i»0,...n-l,  designates  the  error  term  at 
time  t i.  This  error  quantifies  the  data  nonlinearity,  the 
measurement  noise,  and  the  inaccuracy  of  the  derivatives  and 
integrals  of  z\  and  f£.  In  Equation  3-9  jz^|  and  (f 2>  will  be 
treated  as  known  quantities  which  can  be  measured  during  an 
experiment.  The  error  vector  in  Equation  3-9  thus  becomes  a  func¬ 
tion  of  the  system  parameters,  la). 

The  next  step  is  to  Identify  the  parameters  of  the  system 
model.  A  least  squares  approach  is  adopted  In  this 
investigation. 

An  overall  measure  of  the  mismatch  between  the  measured  data 
and  the  model  In  Equation  3-8  Is  established  as  follows. 

e2  «  U)T  {e}  *  ({a)T  fz/l1  -  (f  ,)T  Vfz/i  <a>  -  (f, )\ 
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This  is  referred  to  as  the  squared  error  between  the  measured 
data  and  the  system  model.  This  error  can  be  minimized  through 
the  proper  choice  of  the  parameter  vector  {ah  This  can  be  done 
by  letting 


it2  3e2  it2 

— *r  u  "7  •  •  •  b  9 

3ao  3a  i  3aM+l 


(3-11) 


and  solving  this  sequence  of  equations  for  the  aj,  j*0,...M+l. 
This  solution  can  be  executed  and  the  result,  in  vector  form,  is 


<*>  *  ([V]T  [Zf])'1  [Zf]T  <V  (Z-12) 

The  parameter  vector  chosen  above  is  the  best  estimate  of 
the  system  parameters  in  a  least  squares  sense.  If  the  quantity, 
e2/r\  (where  e2  is  computed  using  Equations  3-9  and  3-10),  is 
relatively  small,  then  Equation  3-3  accurately  represents  the 
measured  system,  if  this  is  not  true,  then  the  model  of  Equation 
3-3  is  Inadequate. 

The  e2  will  be  equal  to  zero  if  the  measured  data  are 
noise  free,  the  system  is  linear,  and  the  computed  derivatives 
and  integrals  of  the  measured  data  are  exact.  Failure  to  meet 
one  or  more  of  these  requirements  will  cause  z2  to  be  nonzero. 

In  practice,  the  parameter  identification  procedure  outlined 
above  can  only  be  used  effectively  when  there  is  little  or  no 
noise  in  the  measurements.  The  method  is  particularly  effective 
when  the  parameter  M  is  set  to  0  because  this  model  will  not,  in 
general,  require  the  computation  of  derivatives  of  z'  and  f.  How¬ 
ever,  when  noise  is  present  and  M  *  1,  the  procedure  loses  its 
accuracy.  If  the  measured  raw  data,  z*  and  f*,  $*0,...n-l, 
are  used  to  obtain  the  derivatives  through  numerical  differentia¬ 
tion,  then  the  estimated  values  of 'z*  and  f  may  be  very  poor  due 
to  amplification  of  the  effects  of  noise. 
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In  view  of  this,  alternate  procedures  for  parameter  identi¬ 
fication  In  the  presence  of  noise  must  be  established.  The 
following  chapters  develop  a  frequency  domain  approach  to  the 
identification  of  system  parameters. 


CHAPTER  4 

IDENTIFICATION  OF  PARAMETERS  IN  FREQUENCY  DOMAIN 

4.0  Identification  of  Parameters 

The  problem  of  parameter  identification  can  be  posed  as  one 
class  in  the  broader  topic  of  optimization.  The  object  of  param¬ 
eter  identification  is  to  make  inferences  about  the  real  world 
and  mathematical  models  on  the  basis  of  measured  input  data.  The 
measured  data  in  this  study  were  assumed  available,  and  were 
simulated  to  represent  the  field  data. 

First,  it  is  assumed  that  there  is  no  noise  present  in  the 
measured  data.  Then,  noise  data  are  introduced.  Note  that  the 
measured  response  data  are  given  as  acceleration  values.  This  is 
realistic  since  the  structural  response  acceleration  is  often  the 
measured  quantity  in  an  experimental  test. 

4.1  Second  Order  System 


Now,  consider  the  second-order  model.  Equation  2-4,  Equa 
tion  2-4  can  be  simplified  by  taking 


cj-  31  -  ft 


The  equation  governing  motion  of  the  system  becomes 


mz  +  aiz  +  aoz  *  f 


Fourier  transform  both  sides  to  obtain 


m(iw)^Z(w)  +  ai(iw)Z(u>)  +  aoZ(w)  s  F(u>) 


(4-3 


*  J  2(t)  e  • 


•“  <  U  < 


F(»)  * /  f(t)  e*1' 


dt  -“<«•»<• 


(4-4) 


are  the  Fourier  transforms  of  z(t)  and  f(t),  respectively.  This 
equation  can  be  rearranged  and  combined  with  Z(w)  and  F(u>)  terms 
on  one  side  of  the  equation  to  obtain 


2  F 

(-mar  +  ao  +  aiiw)  =  j 


(4-5) 


Multiply  each  side  of  the  equation  by  its  complex  conjugate  to 
obtain  the  modulus  squared. 


-mu2  +  ao  +  ai  iwl  *  4^-L 

1  M2 


(4-6) 


Evaluate  the  left  side  and  let  |Fj2/|Z|2  equal  Q(<*>)  to  obtain 


2  ‘  o 

(ao  -  nun  )  +  (aiu)  *  Q(«) 


(4-7) 


This  equation  can  be  used  to  identify  the  parameters  of  second- 
order  linear  system.  However,  this  equation  is  exactly  satisfied 
if  and  only  if:  1)  the  system  under  consideration  is  linear;  2) 
all  measurements  are  noise  free;  and,  3)  the  Fourier  transforms 
of  *z*(t)  and  f(t)  used  in  Equation  (4-7)  are  exact.  When  these 
requirements  are  not  satisfied.  Equation  (4-7)  will  include  an 
error  term  (or  noise  term).  This  practical  case  is  usually  the 
one  that  needs  the  most  attention. 

When  noise  is  present  on  the  measured  input  and  response. 
Equation  (4-7)  can  be  written  as 


2  c  2 

(ao  -  mu'  )  +  (aiu)  *  Q(w)  +  e(u) 


(4-8) 


Note  from  Equation  4-2  that  ao  is  the  equivalent  stiffness  and 
ai  is  the  equivalent  damping  for  the  second-order  linear  system. 
Therefore,  ao  is  greater  in  magnitude  than  ai.  In  view  of  this 
and  the  form  of  Equation  4-8,  ao  can  be  estimated  by  noting  the 
frequency  where  Q(“)  +  e(<*»)  is  a  minimum  whenever  the  equivalent 
damping  factor  is  much  less  than  1  (say  less  than  0.2).  This 
will  be  true  in  most  civil  engineering  systems. 

Denote  the  frequency  where  Q(w)  +  e(u>)  is  a  minimum  by 
wm.  Equation  4-8  shows  that,  approximately, 

ao  *  mwm2  (4-9) 

since  the  first  term  on  the  left  is  approximately  zero  when 
q(u>)  +  e(<u)  is  a  minimum.  Substitute  Equation  4-9  into  Equation 
4-8;  this  yields 

m2(u>m2  -  w2)  +  (ai<o)2  *  Q(w)  +  e(u>)  (4-10) 


Now,  it  is  necessary  to  find  the  coefficient  ai  which  minimizes 
the  £(**>). 

The  coefficients  ai  can  be  evaluated  using  a  least-squares 
approach,  where  the  integral  of  over  a  specific  range  of 
frequencies  is  minimized.  Based  on  Equation  4-10,  set 


rb  («u  *  -  * 

“a 


(4-11) 


where  <oa  amd  are  lower  and  upper  bound  frequencies, 
respectively.  This  frequency  band  is  chosen  so  that  the  system 
response  behavior  can  be  fully  characterized.  It  is  anticipated 


that  the  frequency  band  includes  the  natural  frequency  for  a 
linear  or  slightly  non-linear  system.  For  a  highly  nonlinear 
system,  the  characteristic  frequency,  t^,  will  shift.  However, 
the  frequency  band  can  be  located  by  finding  the  frequency  where 
q(w)  +  e(w)  is  a  minimum,  and  selecting  the  frequency  band  around 
this  frequency. 

In  general,  the  lower  frequency,  n>a,  is  located  at  a  point 
where  its  corresponding  Q(<*»a)  +  e(uia)  value  is  about  5  times 
as  great  as  the  minimum  value  of  Q(id)  +  e(u>);  and  is  the 
higher  frequency  where  Q(^b)  +  e(“b)  is  about  5  times  greater 
than  the  minimum  value  of  Q(w)  +  e{w).  This  method  is  used  to 
select  wa  and  based  on  an  approximate  linear  analysis. 

When  «a  and  are  chosen  in  this  manner,  the  interval  (u>a, 
u>b)  will  be  approximately  the  half  power  bandwidth  of  the  sys¬ 
tem.  This  frequency  interval  reflects  the  characteristics  of  the 
system. 

Take  the  first  partial  derivative  of  e2  with  respect  to 
ai2,  and  set  it  equal  to  0.  The  result  is 


3e2 

H? 


+  m 


Q(w)) 


2 

u  do) 


0  (4-12) 


Simplify  this  to  obtain 


% 

f  ^Q(“)  -  m2 


( 4-12  a) 


Integrate  the  equation  where  possible  to  get 


Let  “a  =  9a  % 


in,  =  g,  u 
b  Hb  m 


(4-14 


where  qa  is  a  coefficient  less  than  one  and  qb  is  a 
coefficient  greater  than  one.  This  equation  can  be  further 
simplified 


This  equation  provides  the  value  of  aj.  'This  is  the  best  esti¬ 
mator  in  the  least-squares  sense.  All  the  parameters  in  this 
equation  are  known  except  the  integral  of  the  w2Q(u>)  term  which 
can  be  evaluated  numerically. 


4.2  Third  Order  Equation 

The  second-order  linear  ordinary  differential  equation  may 
not  be  considered  an  accurate  representation  of  the  hysteretic 
system.  It  is  hoped  that  the  third-order  linear  system  may 
improve  the  accuracy  in  some  sense. 


Equation  2-6  can  be  simplified  by  taking 


Then  Equation  2-6  becomes 


ma2  T  +  mz  +  aiz  +  aoz  =  f  +  a2f 


Fourier  transform  both  sides  to  get 


3  2 

ma2  (i“>)  +  m(iw)  +  ai  (iw)  +  ao 


/ 1  +  a2  (i“ 


The  symbols  used  in  this  equation  have  the  same  meaning  as  in 
earlier  equations.  Multiply  each  side  of  the  equation  by  its 
complex  conjugate  to  yield  the  modulus  squared 


3  2 

ima2“  -  mu>  +  iaiw  +  a2  ipi2 
1  +  ia2  “  ^  | Z | 


(4-19) 


Evaluate  the  left  hand  side  and  let  |F|2/|Z|2  be  replaced  by 
Q(w)  to  obtain 


2  c  2  2  2 

(ao  -  mu>  )  +  u  (ai  -  mw  a2) 

1  +  (a2“j)^ 


=  Q  (u) 


(4-20) 


This  equation  governs  a  third-order  linear  system  in  the  fre¬ 
quency  domain.  Measured  data  can  satisfy  this  equation  exactly 
if  and  only  if:  1)  the  system  under  consideration  is  linear;  2) 
all  measurements  are  noise  free;  and,  3)  the  Fourier  transforms 
used  to  define  Q(w)  are  exact.  These  conditions,  however,  are 
not  usually  met.  In  fact,  the  purpose  of  this  investigation  is 
to  use  the  higher-order  linear  system  to  represent  an  hysteretic 
system.  Therefore,  measured  data  do  not  usually  satisfy  the 


above  equation.  To  account  for  this  explicitly.  Equation  4-20  is 
written 


2  2 
2  c  2  2  c 

(ao  -  mu  )  +  u  (ai  -  mura2) 

- o -  =  Q(<*>)  +  e(w)  (4-21) 

1  +  (a2u) 

e(u)  is  a  noise  term  which  must  be  minimized  by  the  proper  choice 
of  system  parameters.  This  equation  can  be  used  to  identify  the 
system  parameters  following  several  approaches.  Two  of  these  are 
summarized  below.  One  approach  approximates  certain  terms  in 
Equation  4-20  to  obtain  estimates  for  the  system  parameters, 
while  the  other  approach  uses  a  search  technique  to  estimate 
system  parameters. 

The  first  method  to  be  investigated  is  an  approximate  tech¬ 
nique.  When  a2  is  small  compared  to  the  characteristic  fre¬ 
quency,  ^  of  an  SDF  system,  the  (a2u)2  term  in  the  denomi¬ 
nator  can  be  neglected.  This  is  usually  true  when  nonlinear 
deformation  is  hot  too  large.  Eliminate  the  (a2“)2  term  in 
Equation  4-21  to  obtain 

2  2 

(ao  -  mu2)  +  u2  (ai  -  mu2a2)  *  Q(u)  +  e(w)  (4-22) 

When  ai  and  a2  are  small  compared  to  ao  (which  is  usually  the 
case,  the  minimum  of  the  left-hand  side  occurs  near  the  frequency 


(4-23) 


As  previously  described,  the  characteristic  frequency,  can 
be  found  when  the  Q(w)  +  e(u)  term  is  a  minimum.  In  terms  of 
u,,,,  ao  can  be  written 


Expand  the  second  term  on  the  left  side  of  Equation  4-22  and  use 
the  result  of  Equation  4-24.  Neglect  the  m2w^a22  term.  Then 
Equation  4-22  becomes 


m2  (“m2  -  \2)  +  “>2  (a i2  -  2mu>2aia2)  ■  Q(«)  +  e(«)  (4-25) 


Let  ai2  =  bi  and  aia2  =  b2,  then 


(4-26) 


“2  (b i  -  2mw2b2)  -  .^Q(“)  -  m2(wm2  -  w2)^=  e(u)  (4-27) 

When  this  expression  is  evaluated  at  the  discrete  frequency 
w  =  w^,  the  result  is 


“k?  <bl  -  2m“k2b2)  '  (»k  -  ^'“m2  '  “k2)2)  ■ 


*k  (4-28) 


The  system  parameters,  ai  and  a2,  can  now  be  identified  as 
those  which  minimize  the  sum  of  the  squares  of  the  terms  in 
Equation  4-28.  Consider  a  sequence  of  discrete  frequencies 
uniformly  spaced  in  the  interval  (wa,  wfa).  These  are  «*>a, 
o)a  +  aw,  uia  +  2aw,  etc.  Define 


{b}  =  (bi  b2)T 


-  2mu. 


w  +  aw  -  2m(w  +  aw) 
a  a 


[xi]  =  w  +  2 aw 


2m(w  +  2Aw) 


(4-29a) 


(4-29b) 


(4-29d) 


Note  Aw  is  equal  to  (2*/T)  and  T  is  the  total  duration  of  the 
excitation.  (u>a,  u^)  defines  the  range  of  frequencies  over  which 
the  system  is  analyzed.  As  in  the  identification  of  parameters 
of  the  second-order  linear  system,  only  a  portion  of  the  fre¬ 
quency  range  is  used  in  the  parameter  identification. 

This  frequency  range  can  be  chosen  as  before.  In  terms  of 
the  matrices  defined  above.  Equation  4-27  can  be  written  at 
discrete  frequencies  as 

[Xj]  {b}  -  {x2}  *  {e}  (4-30) 

The  vector  {b }  can  be  found  by  minimizing  the  sum  of  the  squares 
of  {e}.  This  is  e2  =  (c)^  {e}.  The  vector  {b }  which  mini¬ 
mizes  e2  is 

(b)  -  (Cx13tC*13»“1C*13t  f*2)  (4-31) 


When  bi  and  b2  have  been  computed,  ai  and  a2  can  be  found  from 
Equation  4-26. 


ai  *  /5T,  a2  =  b2/ai 


(4-32) 


This  solution  provides  approximate  values  for  the  parameters  in 
the  third-order  linear  system  which  represents  the  hysteretic 
system. 
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Another  approach  can  be  used  for  the  estimation  of  param¬ 
eters  in  the  higher  order  linear  system.  This  is  a  search  pro¬ 
cedure  which  iteratively  estimates  the  parameter  values.  Equa¬ 
tion  4-21  can  be  rewritten 


2  2 
(ao  -  mu  )  +  w  (ai  -  mu  &2) 

e(u)  =  Q(u)  -  - 2 -  (4-33) 

1  +  (32“) 

This  quantity  is  a  measure  of  the  mismatch  between  the  measured 
data,  reflected  in  Q(u),  and  the  model,  reflected  in  the  second 
term  on  the  right  side  of  Equation  4-33. 

This  mismatch  can  be  either  positive  or  negative  and  can  be 
used  to  define  one  measure  of  the  difference  between  the  model 
and  the  measured  data  over  a  range  of  frequencies.  This  measure 
is 


e2  *  2  e(u  )  (4-34) 

k  K 

where  the  sum  is  taken  over  those  discrete  frequencies  in  the 
interval  (ua,  ujj).  This  is  the  square  error  of  the  model. 

This  error  is  minimized,  however,  when  the  model  parameters  are 
chosen  to  satisfy  the  sequence  of  equations 


3e2  0  3e2  _  3e2 
dao  u  "Tai  Sa2 


(4-35) 


The  parameters  ao,  ai,  and  a2,  which  satisfy  these  equations 
establish  a  model  which  is  optimal  in  a  least  squares  sense. 

Equations  4-35  can  be  solved  numerically  using  a  search 
technique.  A  computer  program  was  written  to  solve  Equation 


4-35.  The  program,  included  in  the  appendix,  uses  Newton's  meth 
od  to  search  for  the  solution.  The  analysis  procedure  followed 
in  the  computer  program  is  identical  to  that  used  in  solution  of 
the  problem  summarized  in  the  following  section.  The  steps  in 
the  solution  procedure  are  listed  at  the  end  of  Chapter  5. 


CHAPTER  5 


SYSTEM  WITH  TIME-VARYING  PARAMETERS 

5.0  Time- Varying  Parameters  Model 

Sructures  may  exhibit  time-variant  nonlinear  response  to 
strong  motion  excitation.  Time-varying  structural  properties 
were  not  considered  in  the  previous  chapters.  In  this  chapter,  a 
structure  is  modeled  as  a  time  variant  si ngle-degree-of -freedom 
(SDF)  oscillator,  and  a  methodology  is  introduced  to  determine 
its  parameters  using  the  observed  data.  It  is  important  to 
introduce  a  technique  which  can  be  applied  when  noise  is  present 
in  the  measured  data. 

To  demonstrate  this  procedure,  consider  an  SDF  linear  system 
with  mass  m.  Let  the  damping  and  stiffness  parameters  for  this 
system  be  time  varying.  Its  equation  of  motion  is 

mz  +  C(t)z  +  K ( t )  z  *  f  (5-1) 

in  which  z  is  the  displacement  response  of  the  system;  C(t)  and 
K(t)  are  time  variant  damping  and  stiffness  functions  of  the  sys¬ 
tem,  respectively;  and  f  is  the  forcing  function.  It  is  proposed 
that  this  equation  be  used  to  model  the  behavior  of  a  system 
governed  by  Equation  2-1.  Observe  the  system  from  time  0  to  T 

and  assume  that  z{0)  =  0  and  z(0)  3  0. 

The  functions  (C(t)  and  K(t)  are  assumed  to  have  the  form 

C(t)  *  (1  +  at)  co 

K(t)  =  (1  +  St)  ko  (5-2) 

Here,  co  and  ko  are  the  damping  constant  and  the  stiffness  con¬ 
stant,  respectively,  a  and  B  are  constant  coefficients  which  are 
usually  much  less  than  one. 
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In  many  practical  cases,  it  is  observed  that  the  structure 
displays  an  increase  in  damping  and  a  decrease  in  stiffness  when 
the  structure  excites  an  inelastic  response.  This  implies  a  is  a 
positive  constant  and  8  is  a  negative  constant. 

In  this  study,  although  a  and  6  will  be  considered  as  small 
values,  they  will  be  large  enough  to  influence  the  system's  prop¬ 
erties.  This  permits  treatment  of  Equation  5-1  as  a  perturbed 
differential  equation.  When  a  and  8  are  both  equal  to  0  (unper¬ 
turbed),  the  equation  5-1  is  simply  a  second-order  differential 
equation  with  constant  coefficients  which  can  be  easily  solved. 

The  solution  of  Equation  5-1  can  be  written  in  the  form  (for 
example  reference  82) 

2  *  zo  +  a  za  +  8  zB  +  high-order  terms  (5-3) 

Since  a  and  8  are  small,  the  high-order  terms  will  be  neglected. 
Substituting  Equation  5-3  into  5-1  and  expanding  yields 

m(20  +  “Za  +  8?g)  +  CO  (1  +  at)  (20  +  aza  +  Bzg) 

+  ko  (1  +  St)  (20  +  aza  +  B2g)  *  f  (5-4) 

Moving  the  force  term  to  the  left  side  of  the  equation,  grouping 
coefficients  of  the  terms,  1,  a,  and  8,  then  equating  the  coeffi¬ 
cients  to  2ero  results  in 

mzo  +  cozo  +  kozo  *  f  (5-5a) 

mz0  +  coza  +  koza  =  -  cotio  (5-5b) 

mi’e  +  coz8  +  kozB  *  -  kotzo  (5-5c) 

These  equations  approximately  govern  the  system's  response  when 

time  variation  of  the  parameters  is  linear,  as  shown  in  Equation 


5-2.  If  the  excitation  and  the  system's  response  are  known,  then 
these  equations  can  be  used  with  a  time  domain  parameter 
identification  procedure  to  estimate  the  system  parameters. 
However,  when  a  time-domain  parameter  identification  approach  is 
used,  problems  arise  if  noise  is  present  in  the  measured  input 
and  response  signals  (see  Reference  8). 

A  frequency  domain  approach  to  the  identification  of  system 
parameters  is  pursued.  Therefore,  the  equations  of  motion  are 
transformed  to  the  frequency  domain. 

Let 


ZoH 

zo(t)  e~'wt  dt 

(5-6a) 

■/  za(t)  e'1*  dt 
.00 

(5-6b) 

• 

■/  Zg(t)  e"1Ut  dt 

(5-6c) 

define  the 

Fourier  transforms 

of  Z0(t),  za(t),  and  zg(t). 

And  let 

z0(t) 

-  £  /  *.(-)  ‘iut 

du 

( 5— 7a) 

2a(l) 

*  £  l ».(-)  *lBt 

du 

(5- 7b) 

zp(t)  a  7i/z3(“>  e1Ut  d“ 


( 5-7  c ) 


define  the  inverse  Fourier  transforms.  Then  the  Fourier  trans¬ 
form  of  Equation  5-3  is 


Z(«)  =  Zo  +  <*Za  +  6Z&  (5-8) 

It  can  be  shown  that  the  Fourier  transform  of  Equations  5-5a 
through  5-5c  are  given  by 

-mu  Zo(“)  +  co  i<*>Zo(“)  +  ko  Zo{“)  *  F(w)  (5-9a) 

-mu^Za(u)  +  co  i“Za(u)  +  koZa(u)  =  co/zo(w)  +  uZ'o(<ji,)j 


(5-9b) 

-mw2ZB(a>)  +  co  i“Zp(u)  +  koZg(w)  =  -i  ko  Z'  o(**»)  (5-9c) 

Now  solve  the  Equations  5-9  simultaneously.  The  result  is 

Zo(«)  ■  H(u>)  F(u)  (5-10a) 

Za(«)  *  coH(u)  ^H(«)F(»)  +  u>(H'(<o)F(o»)  +  H(u)F'(u))j  (5- 10b) 
Zg(w)  -  -ikoH(u)  (h'F(ui)  +  H(«)  F '(<*>))  (5-10c) 


where  H(u)  is  frequency  ^response  function  and  H'(w)  is  its  first 
derivative.  These  can  be  written  in  the  forms 

2  -1 

H(w)  =  [(ko  -  mur)  +  i(wco)]  (5-lla) 

2  “2 

H'(u)  =  (2mu>  -  ic o)  [(ko  -  mu;  +  i(uco)]  (5-llb) 

F(u>)  is  the  Fourier  transform  of  f ( t)  and  F'(u)  is  its 
derivative. 


F(«)  »/  f(t) 


e-1wt  dt 


(5- 12a) 


dt 


(5— 12b) 


F'  (<u)  =  -i  J  tf  (t)  e_1“ 

Substitute  the  results  from  Equations  5-10  into  Equation 
5-8.  This  yields 

Z(U>)  =  H(w)F(w)  +  oco(h(u))  H(u>)F(u>)  +  uj( H *  (u>)F( o>) 

+  H(u)F'  ( w))+  (-ik  03)  H( w)(h‘  (uj)F( w)  +  H(u>)F' (<*.))  (5-13) 

This  is  the  approximate  frequency  domain  expression  for  the  solu¬ 
tion  of  Equation  5-1.  It  is  considered  accurate  when  both  a  and 
3  are  small.  The  displacement  response  also  can  be  obtained  by 
inverse  Fourier  transformation  of  the  Equation  5-13.  Equation 
5-13  is  used  in  the  identification  process.  Its  use  will  finally 
lead  to  the  estimation  of  the  parameters  from  a  sequence  of  mea¬ 
sured  data. 

5.1  Identification  Procedure 

The  method  described  above  provides  the  solution  for  the 
Equation  5-1  in  the  frequency  domain.  When  the  measured  values 
of  f(t)  are  used  to  estimate  F(w)  and  the  result  is  used  in  Equa¬ 
tion  5-13  to  obtain  Z(u),  this  Z(«)  will  not,  in  general,  match 
the  Z(w)  estimated  from  the  measured  Z(t).  Moreover,  the 
calculated  |Z(w)|  will  not  match  the | Z ( <*>) j  obtained  from  measure¬ 
ment.  The  reasons  for  this  mismatch  are  that  (1)  noise  is  inevi¬ 
tably  present  in  the  measured  input  and  response.  (2)  the  mathe¬ 
matical  model  is  linear,  yet  the  measured  data  come  from  non¬ 
linear  structures  and  (3)  the  discrete  Fourier  transform  of  a 
time  series  is  used  to  represent  the  continuous  Fourier  trans¬ 
form.  In  the  following,  a  brief  theoretical  background  is  pre¬ 
sented  together  with  the  simple  description  of  the  procedure  for 
finding  the  unknown  parameters. 


An  equation  defining  the  mismatch  between  the  measured  data 
and  the  model  of  Equation  5-1  can  be  established.  Let 
|2(m) (u>) |  be  the  modulus  of  the  Fourier  transform  of  the 
measured  structural  response  data.  Let  |Z(w)|  be  the  modulus  of 
the  function  obtained  when  the  Fourier  transform  of  the  measured 
input  data  is  used  in  Equation  5-13.  The  difference  between 
these  function  is  defined 

*(»)  «  |Z(«)|  -  |ZM(m)|  (5-14) 

When  the  discrete  Fourier  transform  is  used  to  approximate 
the  continuous  Fourier  transform  of  a  measured  or  theoretical 
signal,  it  is  defined  at  a  discrete  set  of  frequencies, 
a)k  *  kAo>,  k  =  0,1,..., n-1.  Here  n  and  Aw  relate  to  the  time 
signal  z(t)  and  its  discretization.  It  is  assumed  that  z(t)  is 
available  on  the  interval  (0,  T)  and  is  represented  by  the 
discrete  set  of  values  Zj,  j*0,.../n-l.  Thus,  n  is  the  number 
of  points  vrfjere  the  signal  is  represented.  Au>  is  given  by  2*/T. 
At  a  particular  frequency  <*>  *  «k,  Equation  5-14  becomes 

Ek=  |Zk(»)|  -  |Zk(n,)M|  (5-15) 

ek  can  be  positive  or  negative.  A  quantity  which  is 
always  non-negative  and  which  summarizes  the  differences  between 
the  measured,  jZ^  |,  and  the  theoretical,  jZk|,  structural 
responses  in  the  frequency  domain  over  a  range  of  frequencies  is 
given  by 

e2  =2  e  2  (5-16) 

k  K 

This  is  the  square  error  between  measured  data  and  the 
model.  This  error  can  be  minimized  by  properly  choosing  the 
parameters,  ko,  co,  «,  and  6.  A  method  that  chooses  the 
parameters  this  way  is  a  least  squares  method. 


The  range  of  index  values,  k,  over  which  the  above  sum  is 
taken,  is  not  specified  in  Equation  5-16.  Equation  5-16  need  not 
be  summed  from  0  to  n.  Rather,  the  summation  should  be  carried 
out  over  the  range  of  frequencies  which  includes  those  values  of 
Z|<  containing  significant  information  on  the  behavior  of  the 
system.  In  general,  this  is  the  band  of  frequencies  surrounding 
the  characteristic  frequency  of  the  system. 

Now,  one  can  choose  ko,  Co,  <*,  and  B,  as  those  constants 
which  satisfy 


3 e2  _  3e2  _  3e2  _  3e2  _  n 
dko  3c  o  Sa  3$ 


(5-17) 


The  ko,  co,  «,  and  3,  can  be  located  using  a  search  technique. 

To  simplify  the  analysis,  Newton's  method  is  used  to  minimize 
e2  with  respect  to  ko,  co,  <*,  and  B. 

Newton's  method  converges  very  rapidly  once  an  iterate  is 
fairly  close  to  the  solution.  The  formal  simplicity  and  its 
great  speed  are  the  reasons  why  Newton's  method  is  used  in  this 
study. 

To  assure  convergence  in  the  numerical  analysis,  it  is 
important  to  choose  the  initial  iterate  properly.  A  more 
detailed  discussion  of  the  numerical  procedures  will  be  presented 
later,  in  the  numerical  examples. 

The  steps  in  the  numerical  analysis  are  as  follows: 

1.  Make  the  initial  guesses  at  the  parameter  values,  ko, 

Co,  <*,  and  B. 

2.  Choose  the  computation  increments  Ako,  Aco,  Aa,  and  A3. 


3.  Choose  the  desired  accuracy  measure  (used  to  judge 
convergence) . 
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4.  Compute  the  partial  first  and  second  derivatives  of  tZ 
with  respect  to  ko  using  central  difference  formulas. 

5.  Use  Newton's  method  to  minimize  with  respect  to 
ko. 

6.  Repeat  steps  4  and  5,  this  time  minimizing  with  respect 
to  co,  then  a,  then  £. 

7.  Check  the  result  to  convergence. 

a.  If  convergence  has  occurred,  then  stop  the 
analysis. 

b.  If  convergence  has  not  occurred,  then  repeat  steps  4 
through  6. 

A  computer  program  to  execute  the  procedure  described  above 
has  been  written.  This  program  is  named  PUR  and  a  listing  is 
included  in  the  Appendix. 


CHAPTER  6 

IDENTIFICATION  OF  MULTI-DEGREE-OF-FREEDOM  SYSTEM 


6.0  M.D.F.  System  Model 

The  restoring  force  model  developed  in  the  Equation  2-2  of 
Chapter  2  can  be  extended  for  use  in  the  modeling  of  a  multi- 
degree-of-freedom  (MDF)  system.  In  the  present  investigation  a 
shear  beam  type  MDF  structure  will  be  considered.  Figure  6-1 
shows  the  type  of  structure  under  investigation.  In  this  system 
the  mass,  m^,  is  connected  to  masses  m^  and  by  the 
elements  denoted  k^  and  kj+r 

The  equation  of  motion  governing  the  response  of  a  linear 
mass  excited  undamped  MDF  system  is 


[m]  Vi)  +  [k]  (z)  =  (f } 


(6-1 


where  m  and  k  are  the  mass  and  stiffness  matrices,  respecti vely, 
and  are  written  in  the  form 


n  is  the  number  of  degrees  of  freedom  of  the  MDF  system  in 
in  Figure  6-1.  The  restoring  force  in  this  expression  is  [k]  {z} 
and  this  expression  is  valid  as  long  as  the  response  remains 
linear.  When  the  response  is  inelastic  then  [k]  (z)  must  be 
replaced  by  a  vector  reflecting  the  time  dependent  characteristic 
of  the  restoring  force. 

In  the  present  application  it  will  be  convenient  to  rewrite 
the  equation  of  motion  in  terms  of  relative  displacements.  Let 
yj,  j=2,...n,  denote  the  relative  displacement  between  the 
j-lth  and  jth^  masses,  and  let  yj  denote  the  relative  displace¬ 
ment  between  mass,  mj,  and  the  ground.  Then  the  relation  can 
be  written 


(z)  =  [A]  ly) 


(6-2) 


where 


[A]  = 


10  0 
1  1  0 
1  1  1 


(6-3) 


L i  i  1  •  •  •  ij 

Equation  6-2  provides  an  expression  for  iz),  and  this  result  can 
be  used  in  Equation  6-1  to  obtain  an  alternate  expression  for  the 
equation  of  motion.  This  is 


[M]  Cy)  +  [K]  (y)  *  (f) 


(6-4) 
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where 


CM]  «  Cm]  [A] 


mj  0  0  * 

m£  ^2  0  * 
m3  m3  m3 


( 6— 4a) 


mn  mn  mn 


kj  -k2  0 

0  k2  -k3  0 

0  0  k3  -k4 


CK]  «  [k]  [A]  -  • 


0  0  0 


(6-4b) 


The  term  [k]  {y}  in  Equation  6-4  can  be  expanded  to  obtain  the 


expression 


[K]  = 


klyl  -  ^2 
k2y2  "  k3y3 
k3y3  “  k4y4 


(6-5) 


Note  that  the  restoring  force  applied  to  mass  j  is  k jYj-k .+1, 
The  objective  in  this  investigation  is  to  replace  the  simple  J 


linear  restoring  force  with  a  high  order  linear  model.  Therefore 
in  Equation  6-5  each  of  the  terms  will  be  replaced  with  a 
term  u^.  Each  term  u.  is  governed  by  a  differential  equation. 

J  J 

Let  [U]  and  TR)  be  defined  as  follows 


Then  the  differential  equation  governing  the  motion  becomes 

[M]  (y)  +  [U]  {R>  =  {fJ  (6-8) 

Each  u^,  i*l,...n,  is  governed  by  an  equation  of  the  form 
M 

^  '  C*.M+lV*1  <6'9> 

j-0 

where  qj,  isl,...n,  j=0,...,M+l,  are  the  parameters  which 
characterize  the  restoring  force.  The  superscript  (j)  refers  to 
the  jth  derivative  of  u-j  with  respect  to  time.  This  equation 
is  analogous  to  Equation  2-2  for  the  SDF  system. 
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Equation  6-9  can  be  written  in  matrix  form  using  the 
following  expression.  Let 


Then  the  sequence  of  Equation  6-9  can  be  expressed 

(6-11) 


(6-12) 


To  obtain  the  differential  equation  approximately  governing 
the  motion  of  an  MDF  system,  it  is  necessary  to  combine  Equation 
6-8  and  6-11.  Rearrange  Equation  6-8  to  obtain 


£  [Cj]  ju(j)j  =  (y)  ♦  <y) 

where  {u(j)}  is  the  jth^  derivative  of  the  vector 


[U]  (R )  ■  (f )  -  [M]  (y) 


(6-13) 
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This  is  the  equation  assumed  to  govern  motion  of  the  inelastic 
MDF  system.  The  element  restoring  forces  of  the  MDF  system  are 
governed  by  Equation  6-9.  The  governing  equation  for  the  MDF 
system  is  analogous  to  Equation  2-4,  the  governing  equation  for 
the  SDF  system. 

When  the  parameters  appearing  in  Equation  6-19  are  known  and 
the  input  is  specified,  the  equation  of  motion  can  be  solved  and 
its  solution  provides  the  structural  response  approximation. 
Several  measurements  of  structural  responses  can  be  evaluated. 

For  example,  the  response  displacement,  velocity  and  acceleration 
at  all  structural  degrees  of  freedom  can  be  formed.  Beyond  this, 
other  measures  of  structure  response,  such  as  energy  dissipated 
by  the  structure,  can  be  determined.  Later,  numerical  examples 
of  structural  response  computation  will  be  given. 

The  objective  of  this  investigation  is  to  perform  structural 
parameter  identification  for  inelastic  structures.  Equation  6-19 
could  be  used  to  execute  such  an  identification,  using  measured 
data,  but  the  presence  of  the  higher  derivatives  of  both  input 
and  response  precludes  the  practical  application  of  this 
approach.  (See  Reference  8).  A  better  approach  to  the  parameter 
identification  problem  is  in  the  frequency  domain.  Therefore, 
Equation  6-19  is  transformed.  Fourier  transformation  of  both 
sides  of  Equation  6-19  yields  the  expression 


(  2 IX,]  C°]  M  +  *  [i] 

V=0  / 


tv) 


■  M 

JCCj][D]  (1“)j 

j=0 


(F  > 


(6-20) 


{ Y(«t>) }  =  /  (y(t) )  e'1Uitdt 


(6-20 


(F(w) }  =/*  (f(t)}  eiut  dt 


(6-20 


To  use  this  equation  for  parameter  identification,  the 
following  notation  is  established.  Let 


[H<-)]  *  [D]  tM]  (i»)j+2^  (1»)  CC^j]  +  [I] 


Z  CC,3  [D]  (1»)J 

j=0  J 


(6-21 


The  components  in  [H(w)]  can  be  denoted 


[H(u>)]  = 


Hll(w)  H12(u) 
H2i(“)  H22(u>) 


Hln  ( “) 
«2n  <  <“) 


Hnl<">  Hn2<“>  ’  •  1  Hnn<“> 


(6-2 


Then  Equation  6-23  can  be  expressed  as 


(Y)  =  |H | (F} 


(6- 


The  jth  term  in  the  vector  (Y }  is 


n 

Yj(“)  Hjk(“)Fk(a,)  (6-24) 

k=l 

This  is  the  Fourier  Transform  of  the  structure  response  at  the 
jth  degree  of  freedom  of  the  structure.  The  modulus  of  the 
response  at  this  degree  of  freedom  is 

n 

IVu)l  =  I  2  Hjk(u)  Fk(ui)l’  —  (6-25) 

k=l 


6.1  Identification  Procedure 

Equation  6-25  can  be  used  to  execute  the  identification  of 
the  structural  model  parameters.  If  the  structural  system  under 
consideration  is  linear  with  governing  equation  given  by  Equation 
6-19,  and  if  its  input  and  response  can  be  measured  exactly, 
without  noise,  and  if  the  Fourier  Transform  of  these  signals  is 
performed  exactly,  then  measured  data  can  satisfy  Equation  6-25. 
But,  in  general,  measured  data  are  noisy.  The  Fourier  Transform 
used  in  practical  computation  is  the  Fast  Fourier  Transform 
(FFT) .  And  the  model  of  Equation  6-19  does  not  precisely  charac¬ 
terize  the  structural  system.  Therefore,  Equation  6-25  will  not 
generally  be  exactly  satisfied  where  Yj(<*>)  and  F|c(u>)  are 
obtained  using  measured  data.  In  view  of  this  the  following 


expression  can  be  written 

n 

CJ<“>  ■|V“>I  -  |S  Hjk(»)FkH| 


This  is  an  error  term  and  it  characterizes  the  differences 
between  the  measured  response  data,  and  the  response  that  would 
be  predicted  Ly  the  model.  Equation  6-20. 
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A  least  squares  approach  can  be  used  to  minimize  this 
error.  Let 


define  the  squared  error  of  the  model  in  the  frequency  domain  at 
degree  of  freedom  j.  The  integration  is  taken  over  a  range  of 
frequency  values  such  as  wa  through  These  limits  are 
usually  chosen  to  include  the  frequencies  where  the  system  dis¬ 
plays  power. 

A  measure  of  the  model  error  at  all  degrees  of  freedom  is 

2 

obtained  by  summing  the  at  all  points. 


(6-28) 


This  quantity  reflects  the  model  error  in  the  entire  system. 

This  error  can  be  minimized  with  respect  ot  the  model  parameters 
Cij,  i*n,...n,  j=0,...M+l.  by  solving  the  sequence  of  equation 


"ScT.  s  °»  1=1 . n»  j=0 . M+1 

ij 


(6-29) 


for  the  c-jj.  The  solution  yields  the  system  parameters  which 
best  characterize  the  model  in  a  least  squares  sense. 


i 


Because  of  the  complexity  of  the  system  model,  the  computa¬ 
tion  used  to  solve  Equation  6-29  must  be  executed  numerically. 

2 

This  computation  is  a  search  for  the  minimum  value  of  e  in  the 
nX(M+2)  dimension  space  of  system  parameters,  c. .  ,  i=l,...,n, 

1  J 

j=0,...,M.  The  search  can  be  carried  out  in  one  of  two  ways. 
First,  a  search  can  be  executed  wherein  the  quantity  e2  is 
minimized  sequentially  with  respect  to  each  of  the  parameters 
Cij.  For  example,  e2  can  be  minimized  with  respect  to  c^ 
then  Cj2»  then  c^  etc.  The  minimizing  sequence  is  repeated 
as  many  times  as  necessary  to  obtain  convergence  in  tZ. 

The  second  approach  to  minimization  of  e2  -js  the  gradient 
search  technique.  I«  this  technique,  the  gradient  of  e2  is 
computed  at  a  point  in  the  c-jj  space.  This  information  is  used 
to  choose  a  new  set  of  parameters  urtiere  e2  will  be  smaller  than 
its  original  value.  At  the  newly  adjusted  point,  e2  and  its 
gradient  are  recomputed,  and  this  process  is  repeated  until  e2 
is  near  its  minimum. 

In  the  present  investigation  the  former  analytical  approach 
is  adopted.  Numerical  examples  in  Chapter  8  shows  how  a  least 
squares  computer  program  can  be  used  to  identify  the  parameters 
of  an  MDF  structure. 


CHAPTER  7 

ENERGY  DISSIPATED  RELATED  TO  CONCRETE  DAMAGE 
7.0  Introduction 

This  presentation  describes  and  evaluates  an  experimental 
study  of  the  strength  reduction  and  behavior  of  plain  concrete 
subjected  to  cyclic  loading.  It  is  recognized  that  concrete  is 
damaged  by  application  of  stresses  lower  than  the  ultimate 
stress.  The  concrete  fracture  process  begins  at  very  low  stress 
and  is  continuous. 

The  damage  caused  by  loading  to  small  stresses  is  slight  and 
each  subsequent  loading  over  the  same  stress  range  produces  a 
negligible  increase  in  damage.  However,  as  the  loading  stress  is 
increased,  more  damage  occurs.  Stresses  with  peak  values  in  the 
range  of  40  percent  to  100  percent  of  the  ultimate  stress  produce 
considerable  damage  and  subsequent  loading  over  the  same  range 
cannot  be  neglected.  In  practical  situations,  when  a  severe 
excitation  is  applied  to  a  structure,  it  is  not  uncommon  for  the 
peak  stress  to  go  beyond  the  50  percent  level  of  ultimate 
stress. 

When  loading  is  repeated,  damage  accumulates  in  a  concrete 
specimen;  consequently  it  no  longer  retains  its  original 
strength.  This  concept  suggests  that  it  might  be  useful  to 
attempt  a  quantitative  evaluation  of  damage  occurring  in  concrete 
during  the  cyclic  loading.  The  objective  of  this  study  is  to 
demonstrate  that  concrete  damage  and  strength  reduction  are 
related  to  energy  dissipation  under  repeated  loading. 

When  energy  is  dissipated  during  loading  and  unloading,  an 
hysteresis  loop  is  formed  in  the  stress-strain  curve  as  shown  in 
Figure  7-1.  The  area  enclosed  represents  the  total  energy  dissi¬ 
pated  during  one  cycle  of  loading.  This  dissipated  energy  may  be 
classified  into  two  parts,  namely,  damage  and  damping  energy 


dissipation.  However,  the  total  accumulated  energy  dissipation 
is  of  primary  interest  in  this  experimental  study.  A  more 
detailed  discussion  of  the  energy  dissipation  mechanism  is  given 
in  Reference  [69]. 

There  are  many  methods  that  can  be  used  to  detect  and  assess 
damage  in  concrete.  Among  the  most  frequently  used  techniques 
are  those  which  assess  change  in  initial  elastic  modulus,  and 
those  which  measure  acoustic  emissions,  change  in  pulse  velocity, 
and  energy  dissipation.  In  the  present  study,  the  dissipated 
energy  method  will  be  adopted.  Attention  will  be  focused  on  this 
means  for  measuring  damage  because  of  its  intuitive  relationship 
with  the  energy  dissipation  of  an  hysteretic  structure  under 
dynamic  loads.  Other  methods  may  be  considered  as  in  References 
[69,  70,  and  71]. 

In  the  present  investigation  a  sequence  of  physical  experi¬ 
ments  was  performed.  In  each  experiment  a  concrete  cylinder 
(specifications  given  below)  was  loaded  in  uniaxial  compression. 
The  load  applied  to  each  cylinder  was  a  cyclic  load,  and  the 
energy  dissipated  was  calculated.  This  was  done  by  plotting  the 
stress  versus  strain  diagram  and  by  determining  the  area  enclosed 
within  the  hysteresis  loops.  Varying  amounts  of  energy  were 
dissipated  in  the  various  test  specimens,  and  upon  completion  of 
cyclic  testing  each  specimen  was  loaded  to  failure  in  order  to 
determine  its  residual  strength.  For  each  test  specimen,  dissi¬ 
pated  energy  and  residual  strength  were  recorded  and  the  relation 
between  these  quantities  was  established. 

More  details  of  the  testing  procedure  are  given  below, 
together  with  fresh  concrete  properties  observed  during  the  mix¬ 
ing.  These  may  provide  a  useful  reference  for  the  concretes  used 
in  the  test. 
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The  concrete  specimens  tested  in  this  investigation  have  the 
mix  details  and  plastic  properties  of  fresh  concrete  as  shown  in 
Table  7-1. 

TABLE  7-1. 

Concrete  Mix  Details 

Aggregate  Ratio  of  Maximum 
Cement  Coarse:Fine  Aggregate 
Type  of  Cement  W/C  Ratio  Ratio  Aggregate  Size 

(in) 

Type  1  A  0.53  4.8  60:40  3/4 

Plastic  Properties  of  Fresh  Concrete 


Mix  No. 

Slump 

(in) 

Air 

% 

Temperature 
(degrees  C) 

Weight 

lb/ft* 

2  WE 

4 

3.5 

27 

145.8 

3  WE 

4  1/4 

4 

30 

148.96 

The  specimens 

were  all 

cast 

in  6- in  x  12- 

in  steel 

molds.  The  concrete  mix  proportions  were  constant  for  all  the 
specimens.  The  specimens  were  tested  at  a  consant  loading  rate 
of  1000  lb/sec  in  a  RIEHLE  compression  testing  machine.  The 
force  versus  strain  results  of  each  test  were  plotted  with  an  x-y 
electronic  recorder.  This  recorder  was  connected  by  an  elec¬ 
tronic  compressometer  which  is  properly  designed  for  this  spe¬ 
cific  purpose  as  shown  in  Figure  7-2.  The  test  machine  was 
properly  calibrated  before  the  test. 


To  ensure  a  uniform  displacement  of  the  specimens,  thin  sul¬ 
fur  caps  on  the  two  end  surfaces  of  the  specimens  were  employed 
and  were  allowed  to  harden  before  testing.  Specimens  were  cured 
in  water  in  the  curing  tank  at  25*C  for  14  days  and  28  days. 


In  assessing  the  energy  dissipated  and  residual  strength, 
specimens  were  subjected  to  a  series  of  cycles  of  loading  and 
unloading.  The  specimens  were  loaded  up  to  a  value  in  the  range 
of  stress,  90  -  94  percent  of  the  ultimate  stress.  This  ensured 
that  damage  occurred  for  every  cycle  of  loading.  At  the  end  of 
each  cylce  of  loading  and  unloading,  the  testing  machine  was 
returned  to  a  rest  position,  and  reloading  was  commenced  immedi¬ 
ately.  To  ensure  that  concrete  characteristics  would  be  as 
nearly  uniform  as  possible,  all  the  tests  in  each  sequence  were 
run  in  one  day. 

7.1  Discussion  of  Results 


Numerous  physical  experiments  were  conducted  in  this  inves¬ 
tigation  and  characteristics  of  concrete  accumulating  damage  can 
be  derived  from  the  individual  tests  and  all  the  tests,  jointly. 
In  the  following  section  the  characteristics  of  individual  tests 
are  discussed  first;  then  damage  characteristics  related  to  the 
entire  test  sequence  are  discussed. 

A  typical  stress-strain  diagram  obtained  during  one  experi¬ 
ment  is  shown  in  Figure  7-3.  A  number  of  characteristic  features 
can  be  extracted  from  this  result.  On  the  initial  cycle  the 
specimen  was  loaded  to  a  stress  near  its  ultimate  (95  to  98  per¬ 
cent).  It  can  be  seen  that  the  most  significant  change  in  behav¬ 
ior  between  consecutive  loading  cycles  occurs  between  the  first 
and  second  cycle. 

The  first  loading  curve  shows  more  curvature  than  the  fol¬ 
lowing  reloading  curves  in  which  curvature  tends  to  diminish. 

The  reloading  curves  show  progressively  decreasing  slopes.  This 
may  be  attributed  to  structural  degradation  of  the  specimen. 

Another  measure  of  degradation  can  be  established  by  plot¬ 
ting  the  initial  elastic  modulus  for  a  particular  cycle  versus 
energy  dissipation  prior  to  that  cycle.  This  is  shown  in  Figure 
7-4.  As  the  energy  dissipated  gradually  increases,  the  initial 
elastic  modulus  diminishes. 


fint  loa^int  md  unkMdmj 
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The  above  discussion  was  based  on  one  typical  sample.  A 
similar  discussion  could  be  given  for  the  other  samples.  Figures 
7-5,  7-6,  and  7-7  show  the  stress-strain  curves  for  some  other 
specimens  tested  during  this  investigation. 

Some  general  characteristics  of  the  accumulation  of  damage 
in  concrete  specimens  can  be  derived  from  the  entire  collection 
of  results.  A  total  of  24  concrete  specimens  were  tested  in  this 
investigation. 

As  noted  earlier,  the  specimens  were  subjected  to  cyclic 
loadings  inducing  different  amounts  of  energy  dissipation  in  the 
various  cylinders.  Not  all  specimens  were  cycled  to  failure.  At 
least,  three  of  the  specimens  were  tested  for  the  determination 
of  the  ultimate  strength.  Other  specimens,  however,  were  cycled 
till  failure.  The  remainder  of  the  specimens  were  cycled  till  a 
certain  amount  of  energy  was  dissipated;  then  these  were  loaded 
to  failure  in  order  to  find  their  residual  strength. 

Using  these  data,  a  characteristic  of  the  specimens  can  be 
extracted.  The  total  energy  dissipated  by  each  particular  speci¬ 
men  was  plotted  against  the  residual  strength  of  the  specimen  as 
shown  in  Figure  7-9.  Another  result  can  be  illustrated  by  plot¬ 
ting  the  total  energy  dissipated  versus  percentage  of  decrease  in 
strength  as  shown  in  Figure  7-8.  Both  diagrams  show  a  decrease 
in  strength  as  the  total  energy  dissipated  is  increased.  Since 
only  a  limited  number  of  specimens  were  tested,  no  direct  mathe¬ 
matical  expression  relating  the  residual  strength  to  the  total 
energy  dissipated  was  obtained.  While  such  a  relation  could  be 
established,  further  testing  is  required  to  derive  a  general 
relationship.  However,  the  present  results  provide  the  informa¬ 
tion  needed  to  conclude  that  energy  dissipation  is  truly  related 
to  the  residual  strength. 


stress-strain  diagram  for  a  concrete  test  specimen 
clic  load.  Peak  stress  90  percent  of  failure 


ypical  stress-strain  diagram  for  a  concrete  test  specimen  under 
:yclic  load.  Peak  stress  90  percent  of  failure  stress  for  Batch  #2 
28  days  curing). 


cal  stress-strain  diagram  for  a  concrete  test  specimen  under 
ic  load.  Peak  stress  90  percent  of  failure  stress  for  Batch  #3 
days  curing). 
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The  experimental  technique  described  above  provides  an 
approach  for  the  estimation  of  damage  in  concrete.  Based  on  the 
physical  experiments,  the  following  conclusions  can  be  made: 

1.  The  most  significant  change  in  the  properties  of  the 
concrete  occur  between  the  first  cycle  and  second  cycle 
when  loading  in  the  first  cycle  is  severe. 

2.  The  initial  elastic  modulus  of  the  specimens  gradually 
diminishes  as  energy  is  dissipated.  This  implies  that 
the  damage  of  concrete  under  cyclic  loading  occurs 
progressively. 

3.  The  energy  dissipated  in  a  concrete  specimen  is 
adversely  related  to  residual  strength.  As  the  energy 
dissipated  increases,  the  residual  strength  decreases. 
Therefore,  energy  dissipated  m^y  be  used  to  predict  the 
damage  of  a  structure  under  a  severe  loading.  Moreover, 
total  energy  dissipated  may  be  considered  as  an 
indicator  of  the  degree  of  damage  in  an  hysteretic 
structure. 

Some  restrictions  apply  to  the  above  conclusions.  The  work 
is  limited  to  behavior  in  uniaxial  compression.  Other  types 
of  loading  are  possible  and  further  tests  are  required  to  charac¬ 
terize  damage  under  general  loading.  It  has  been  assumed  that 
the  creep  effect  is  small  enough  to  be  neglected  in  this 
investigation. 


CHAPTER  8 

NUMERICAL  EXAMPLES 


8.0  Data  Description 

In  this  chapter  numerical  examples  are  presented  which 
demonstrate  the  use  of  the  analytic  procedures  developed  in  the 
previous  chapters.  The  first  set  of  examples  demonstrates  the 
identification  of  the  model  parameters  for  linear  and  hysteretic, 
single-degree-of-freedom  (SDF)  structures  both  when  measurement 
noise  is  and  is  not  present.  One  example  demonstrating  the  time 
domain  approach  to  parameter  identification  is  summarized.  Two 
examples  showing  the  frequency  domain  appraoch  are  presented. 

Another  example  demonstrates  the  application  of  the  analysis 
presented  in  Chapter  6.  The  frequency  domain  approach  is  used  to 
identify  the  parameters  of  a  multi-degree-of-freedom  (MDF) 
structure. 


The  input  used  to  excite  the  SDF  system  in  all  numerical 
examples  is  a  decaying  exponential,  oscillatory  function.  It  is 
generated  using  the  formula 


f (t)  =  e 


•at 


N 

1 

j  =  l 


cos (wjt  -  ♦j) 


0  <  t  <  T 


(8-1) 


where  a,  c.,  j=l,...,N,  and  u>  j=l,...,N  are  constants.  <t>. , 

J  J  J 

j=l,...N,  are  phase  angles  which  are  random  variable  realiza¬ 
tions;  these  random  variables  are  independent  and  uniformly 
distributed  on  the  interval  (0,  2*).  o  is  a  decay  rate.  The 
Cj,  j=0,...N,  are  constants  which  determine  the  amplitudes  of 
the  excitation.  All  the  values  of  cj  are  taken  as  equal  to  c 
in  all  cases  for  the  examples,  wj,  j=l,...N,  are  equally 
spaced  in  the  interval  including  the  characteristic  frequency  of 
the  system  being  analyzed. 


The  forcing  function  defined  above  was  generated  at  discrete 
times.  Specifically,  f(t)  was  evaluated  at  the  times  t=t -|  =1  At , 
1=0,... N-l.  A  computer  program,  named  FORCE,  which  generates  the 
excitation  of  Equation  8-1  was  used  in  these  numerical  examples. 

Three  distinct  signal  types  were  identified  in  the  numerical 
examples.  The  first  type  used  a  computer  program,  named  BILIN, 
to  compute  linear  and  nonlinear  response.  BILIN  can  be  used  to 
find  the  displacement,  velocity,  and  acceleration  response  of  a 
given  bilinear  hysteretic  system  to  an  arbitrary  input.  It  also 
computes  the  energy  dissipated  by  the  structure  during  the 
response.  The  second  type  used  a  computer  program,  named  TIMEVA, 
to  compute  the  response.  This  program  computes  a  linear  time 
dependent  response  defined  by  Equations  5-1  and  5-2  with  a,  P, 
co,  and  ko  constants.  The  third  type  used  a  computer  program 
named  BLNMDF  to  compute  the  MDF  structure  response. 

White  noise  was  used  whenever  measurement  noise  was  added  to 
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the  signals.  The  white  noise  is  normally  distributed,  N(0,on) . 

A  subroutine  named  NOISE  was  used  to  generate  the  noise  signal. 
The  noise  signals  were  added  to  the  generated  input  and  response 
signals  in  the  following  manner.  First,  the  excitation  and 
response  signals  were  generated  using  programs  FORCE,  BILIN, 
BLNMDF,  and  TIMVA.  Then  noise/signal  ratios  were  selected  and 
used  to  obtain  the  variances  of  the  noise  signals.  The  noise 
signals  were  generated  as  sequences  of  independent  random  vari¬ 
ables,  and  directly  added  to  the  excitation  and  response.  These 
noisy  signals  were  then  used  as  inputs  to  do  the  identification. 
Note,  no  filtering  procedure  was  used  on  the  simulated  measured 
signals  during  the  identification  process. 

Three  basic  models  were  used  to  represent  the  hysteretic 
systems.  These  are  the  second  order  linear,  time  invariant, 
third  order  linear,  time  invariant,  and  second  order  linear,  time 


varying  models.  All  the  model  parameters  were  identified.  The 
identification  procedures  and  formulations  were  described 
previously.  Different  identification  approaches  may  be  applied 
for  the  same  model . 

8.1  Example  1 

This  example  solves  a  parameter  identification  problem  using 
the  direct  time  domain  approach,  summarized  in  Equations  3-4 
through  3-12.  The  parameters  of  the  shock  input  are  listed  in 
Table  8.1.  A  typical  forcing  function  history  generated  using 
these  parameters  is  shown  in  Figure  8-la.  The  derivative  of  the 
forcing  function  is  shown  in  Figure  8-lb. 

The  notations  for  the  parameters  used  in  specifying  the 
numerical  examples  are  those  used  in  the  text.  Some  additional 
notations  are  defined  here,  c  is  the  viscous  damping  in  an  SDF 
system,  k  is  the  initial  stiffness  in  a  bilinear  hysteretic 
structure,  ky  is  the  yield  stiffness  of  a  bilinear  h  steretic 
structure.  Zy  is  the  yield  displacement  of  a  bilinear  hyster¬ 
etic  structure.  zmax  is  the  maximum  displacement  of  an  SDF 
structure. 

In  this  numerical  example,  two  basic  problem  types  are 
solved.  These  are  summarized  below. 

(1)  An  input  is  generated  using  Equation  8—1.  The  input  is 
used  to  excite  a  linear  SDF  system  with  viscous  damp¬ 
ing.  The  structural  input  and  response  are  stored,  and 
no  noise  signals  are  added  to  the  input  and  response. 
Then  the  input  and  response  are  used  to  identify  the 
model  parameters,  aj,  j=0,...M+l,  from  Equation  3-5. 

(2)  An  input  is  generated  as  in  1,  above,  but  here  the 
response  of  a  bilinear  hysteretic  system  is  computed. 

No  noise  signals  are  added  to  the  input  and  response. 
The  input  and  response  are  used  to  identify  the  model 
parameters,  a-j,  j=0,...M+l. 


FORCE/ SEC 


[«I«] 


For  these  two  basic  problems  eight  system  identification 
problems  are  solved.  In  cases  one  and  five,  problem  type  one 
(above)  is  considered.  In  case  one,  a  second  order  model  (M=0), 
Equation  2-3)  is  identified.  In  case  five,  a  third  order  model 
(M=l)  is  identified.  In  cases  two  through  four  and  six  through 
eight,  problem  type  two  (above)  is  considered.  In  cases  two 
through  four,  a  second  order  model  (M=0)  is  identified.  In  cases 
six  through  eight,  a  third  order  model  (M=l)  is  identified.  The 
successive  cases  involve  increasing  degrees  of  yielding.  In 
these  eight  cases,  the  parameters  of  the  second  order  model,  ao, 
ai,  and  a2  (M=l)  are  identified. 

Once  the  parameters  of  Equation  3-5  have  been  estimated,  the 
energy  dissipated  by  the  model  is  computed;  and  this  is  compared 
to  the  energy  dissipated  by  the  actual  system  as  computed  in 
BILIN.  This  computation  is  the  one  discussed  in  Chapter  2  and 
given  by  Equation  2-7b.  The  energy  computations  are  performed  in 
programs  ENER2,  for  second-order  systems,  and  ENER3,  for  third- 
order  systems.  The  computations  are  performed  using  an  incre¬ 
mental  form  of  the  governing  Equation  3-3. 

The  responses  of  some  SDF  systems  to  the  shock  input  in 
Figure  8-la  were  computed.  First,  the  response  of  a  linear  sys¬ 
tem  was  computed  for  analysis  in  cases  1  and  5.  The  computed 
displacement  response  is  plotted  versus  time  in  Figure  8-2a.  The 
SDF  structure  spring  restoring  force  versus  displacement  is 
plotted  in  Figure  8-2b.  The  very  slightly  nonlinear  response  of 
an  SDF  structure  was  computed  for  analysis  in  cases  2  and  6,  but 
this  response  is  not  shown.  A  more  nonlinear  response  was 
computed  for  analysis  in  cases  3  and  7.  The  displacement 
response  versus  time  is  plotted  in  Figure  8-3a,  and  the  spring 
restoring  force  versus  displacement  is  shown  in  Figure  8-3b.  The 
first  figure  shows  a  residual  plastic  displacement  as  the 
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Figure  8-2a.  Response  of  a  linear  system  to  the  force  in 
Figure  8 -la.  k  =  ®. 


10000.0 


10000.0  <- 

2.0 


-1.0 


DISPLACEMENT  -  IN 

Figure  8-2b.  Spring  restoring  force  versus  displacement 
for  linear  system. 
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response  vibrations  diminish.  The  second  graph  shows  the  perma¬ 
nent  set  as  lateral  displacement  of  the  horizontal  axis  inter¬ 
cept.  Finally,  a  very  nonlinear  response  was  computed  for 
analysis  in  cases  4  and  8.  The  displacement  response  versus  time 
is  shown  in  Figure  8-4a  and  the  spring  restoring  force  versus 
displacement  is  shown  in  Figure  8-4b.  A  considerable  permanent 
set  is  evident,  as  the  motion  diminishes,  from  the  first  figure. 
The  seond  figure  shows  that  large  plastic  deformations  occur  in 
the  structure  in  both  directions  of  motion. 

The  energy  dissipated  by  each  structural  system  is  listed 
with  the  structural  parameters  in  Table  8.1.  Eg  is  that  energy 
dissipated  due  to  the  action  of  the  inelastic  spring  and  the 
action  of  the  viscous  damper. 

Using  the  forcing  function  input  described  above,  and  the 
computed  responses,  the  parameters  of  the  structural  systems  were 
identified.  The  results  of  the  parameter  identifications  are 
given  in  Table  8.1.  The  energy  dissipated  when  the  identified 
systems  respond  to  the  shock  input  is  listed  in  Table  8.1  next  to 
the  identified  parameters. 

Figures  8-5a  through  8-9b  show  the  computed  responses  of 
some  of  the  identified  systems.  The  figure  titles  indicate  which 
systems  generate  the  responses  shown.  The  top  (or  "a")  figures 
show  the  computed  responses  versus  time.  The  bottom  (or  "b") 
figures  show  the  computed  restoring  forces,  spring  force  plus 
damper  force,  versus  displacements. 

Figure  8-10  compares  three  responses.  These  are:  1)  the 
actual  nonlinear  structural  response  obtained  using  BILIN  in  the 
slightly  nonlinear  cases  2  and  6,  2)  the  response  executed  by  the 
identified  model  described  in  case  2,  and  3)  the  response 
executed  by  the  identified  model  described  in  case  6.  The  two 
latter  responses  practically  overlay  and  very  nearly  equal  the 
actual  response. 


Figure  8-4a.  Response  of  nonlinear  system  to  the  force  in 
Figure  8-la.  k  =  5000 
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Figure  8-4b. 
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Figure  8-7a.  Displacement  resp< 
force  in  Figure  8' 
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Figure  8-8a.  Displacement  response  of  identified  system  to 
force  in  Figure  8-la.  (Case  7) 
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Figure  8-9a.  Displacement  response  of  identified  system  to 
force  in  Figure  8-la.  (Case  8) 
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Figure  8-9b.  Total  restoring  force  versus  displacement  for 
identified  system.  (Case  8) 
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Figure  8-10. 


The  comparison  between  measured  response  (light 
line)  and  identified  response  (dark  line)  for 
both  second-order  (case  3)  and  third-order  (case  7) 
approximate  systems.  (The  identified  responses 
are  practically  the  same.) 


80 


TABLE  8.1.  PARAMETERS  AND  RESULTS  FOR  EXAMPLE  1 
Input  (Structural  Excitation) 


Case  Number 

N 

a 

C 

% 

n 

1  through  8 

50 

i  1.0 

300  32 

.4 

1256.4 

256 

Structure  Parameters 

Case  Number 

m  c 

k 

h. 

zy 

ed 

1.5 

0.259  7.77 

5829 

CO 

23292 

2,6 

0. 

259  7.77 

5829 

0 

7700 

23588 

3,7 

0. 

259  7.77 

5829 

0 

7000 

23674 

4,8 

0. 

259  7.77 

5829 

0 

5000 

21959 

Identification  Parameters  and 

Results 

Case  Number  > 

M 

a0 

al 

a 

2 _ 

ed 

1 

2 

5829 

7.77 

21533 

2 

2 

5825 

7.79 

21533 

3 

2 

4574 

8.13 

21779 

4 

2 

1620 

9.69 

23844 

5 

3 

5829 

7.77 

0.0 

21533 

6 

3 

5823 

7.85 

0.0 

21498 

7 

3 

2527 

99.65 

0.0158 

20824 

8 

3 

610 

108.11 

0.0175 

18666 

Figure  8-lla  compares  the  actual  nonlinear  response  of  cases 
3  and  7  to  the  second-order  (M=0)  model  repsonse  of  case  3. 

Figure  8- lib  compares  the  actual  nonlinear  response  of  cases  3 
and  7  to  the  third-order  (M=l)  model  response  of  case  7.  Both 
models  simulate  response  amplitudes  quite  well;  the  third-order 
model  is  slightly  closer  than  the  second-order  model  in  that  the 
third-order  model  provides  a  slightly  better  phase  match  to  the 
actual  response  than  the  second-order  model. 

In  this  numerical  example,  no  cases  are  included  where  noise 
was  added  to  the  forcing  function  input  and/or  the  acceleration 
response.  Such  cases  were  analyzed,  but  the  results  were  so  poor 
that  they  are  not  summarized  here.  These  results  showed  that  the 
direct,  time  domain  parameter  identification  technique  is  not 
effective  in  parameter  analysis  when  recording  noise  is  present. 

The  numerical  examples  summarized  here  show  that  the  direct, 
time  domain  parameter  identif iction  technique  can  be  used  effec¬ 
tively  when  the  measured  signals  are  noise  free.  This  is  best 
confirmed ‘by  reference  to  Figures  8-10,  8-lla,  and  8-llb.  These 
show  that  the  linear  model  response  can  be  made  to  match  the 
nonlinear  response  well. 

The  energy  dissipation  results,  summarized  in  Table  8.1, 
show  that  the  third-order  models  provide  the  best  simulatior  for 
a  nonlinear  hysteretic  system. 

8.2  Example  2. 

In  this  example,  a  sequence  of  parameter  identification 
problems  is  solved  using  the  frequency  domain  approach.  Three 
programs,  FREQIO,  PUR3,  and  PUR,  were  written  to  execute  the 
parameter  identifications. 

FREQIO  performs  approximate  frequency  domain  parameter 
identification  for  second  and  third  order  linear  models.  It 
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Figure  8-1  la.  The  comparison  between  measured  response 
(light  line)  and  identified  response  for 
second-order  system.  (Case  3) 
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-lib.  The  comparison  between  measured  response 
(light  line)  and  identified  response  for 
third-order  system.  (Case  7) 


accepts  both  an  input  signal  from  FORCE  and  a  response  signal 
from  BILIN  or  TINEVA.  When  desired,  the  white  noise  signals  are 
added  to  the  corresponding  input  data.  Then  FREQID  performs  the 
necessary  Fourier  transforms  and  other  data  operations.  Follow¬ 
ing  this,  the  parameter  identification  is  executed.  One  opera¬ 
tion  required  In  the  parameter  identification  is  estimation  of 
the  characteristic  frequency.  This  can  be  done  simply  by  search¬ 
ing  Q(w)  +  e(u)  for  a  minimum  value.  However,  a  more  precise 
method  for  determining  the  minimum  value  of  Q(«)  +  e(w)  defined 
In  Equations  4-8  and  4-11  involves  use  of  a  least  square  method. 
In  this  improved  method,  Q(w)  +  e(w)  is  used  to  identify  the 
characteristic  frequency  and  other  parameters. 

In  FREQIO,  an  Important  assumption  was  made  for  the  third- 
order  model.  Equation  4-21.  In  particular,  it  was  assumed  that 
a2  is  small  in  value.  The  parameters  of  this  model  may  be  iden¬ 
tified,  without  the  assumption  that  a2  is  small,  by  the  search 
technique.  The  computer  program  PUR3  was  written  for  this  pur¬ 
pose.  A  detailed  description  of  this  method  was  given  In 
Chapter  5. 

Program  PUR  performs  parameter  identification  for  second- 
order  time  varying  linear  systems.  The  approach  is  based  on  the 
procedure  described  In  Chapter  5.  The  program  accepts  the  Inputs 
and  responses  generated  In  the  programs  FORCE  and  BILIN  or 
TIMEVA,  with  or  without  noise.  Four  parameters,  namely  ko,  <*, 
co,  and  6  are  Identified.  The  program  employs  a  search  tech¬ 
nique;  the  Initial  estimators  can  be  chosen  by  using  the  Identi¬ 
fied  parameters  obtained  from  any  of  the  methods. 

Once  the  parameters  have  been  estimated,  the  energy  dissi¬ 
pated  by  the  model  is  computed.  This  result  together  with  the 
predicted  response  is  compared  to  both  the  energy  dissipated  and 
the  response  of  the  actual  system.  The  energy  and  response 
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computations  are  performed  In  program  ENER2  and  ENER3  for  the 
second  and  third  order  system,  respectively. 

From  the  above  description,  the  methods  used  in  the  deter- 
mination  of  the  system  parameters  can  be  summarized  as  follows: 

Method  1.  Performs  parameter  identification  for  the 
second-order  system  in  the  frequency  domain 
utilizing  Equation  4-15.  No  fitting  equation  for 
Q{«)  +  e(u)  is  applied. 

Method  2.  Performs  parameter  identification  for  the  third- 
order  system  in  the  frequency  domain  utilizing 
the  equations  from  4-29a  to  4-31.  No  fitting 
equation  for  Q(«)  +  e(u)  is  applied. 

Method  3.  Performas  parameter  identification  using  the  same 
approach  as  Method  1  except  the  input  data 
Q(u>)  +  e(u)  are  replaced  by  the  fitted  polynomial 
equation.  This  additional  analysis  is  done  in  a 
subroutine  called  FIT. 

Method  4.  Performs  parameter  identification  using  the  same 
approach  as  Method  2  and  using  the  same  procedure 
described  in  Method  3. 

Method  5.  Performs  parameter  identification  for  the  third- 
order  system  using  the  search  method  described  in 
Chapter  5.  The  operation  is  executed  in  a 
program  called  PUR3.  Prior  estimates  obtained 
from  the  above  methods  are  used. 

Method  6.  Performs  parameter  identification  for  the  second- 
order  time-varying  parameter  system  in  the  fre¬ 
quency  domain.  The  search  method  described  in 
Chapter  5  is  used.  This  method  also  requires 
prior  estimators  which  can  be  obtained  from  the 
information  supplied  in  Methods  1  through  5. 


All  the  methods  described  above  can  be  used  to  estimate  the 
parameters  for  the  linear  and  nonlinear  systems  even  when  noise 
Is  present.  The  duration  of  the  excitation  must  be  long  enough 
to  characterize  the  system  parameters. 

In  the  following  numerical  examples,  four  basic  problems  are 
solved.  These  cases  involving  different  degrees  of  nonlinearity 
In  the  system  response  are  summarized  below. 

Case  1.  An  input  excitation  is  generated  using  Equation 
8-1.  The  input  is  used  to  excite  a  linear  SDF 
system  with  viscous  damping.  The  excitation  and 
linear  response  are  used  to  Identify  the  model 
parameters.  Noise  signals  can  be  added  to  the 
generated  input  and  response.  If  required. 

Case  2  An  excitation  Input  Is  generated  as  In  Case  1, 
but  here  the  response  of  a  bilinear  hysteretlc 
system  is  computed.  Yielding  occurs  In  the 
response.  The  degree  of  nonlinearity  was 
designed  using  a  comparison  between  the  yield 
displacement  of  the  bilinear  system  and  the 
maximum  displacement  of  the  linear  system.  Let 
the  yield  displacement  of  the  bilinear  system  be 
Zy.  Let  the  maximum  displacement  of  the  linear 
system  be  z^.  In  this  case,  zmax  is  taken  as 
6.7  and  Zy  Is  equal  to  6.0. 

Case  3.  Same  as  Case  2,  but  Zy  Is  equal  to  5.0. 

Case  4.  Same  as  Case  2,  but  Zy  Is  equal  to  4.0. 

Case  5.  An  Input  excitation  Is  generated  as  in  Case  1. 

The  Input  Is  used  to  excite  a  linear  SDF  system 
with  time  varying  damping  and  stiffness.  The 
excitation  and  response  are  used  to  identify  the 
model  parameters.  Noise  signals  are  added  to  the 
simulated  Input  and  response  when  required. 
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The  example  carries  out  the  parameter  identification  using 
the  methods  described  above.  Specifically,  methods  1  through  6 
are  used  to  identify  the  parameters.  The  parameters  of  the  input 
excitation  are  listed  in  Table  8.2.  The  notation  for  the  param¬ 
eters  was  specified  above. 

TABLE  8.2.  PARAMETERS  OF  THE  FORCING  FUNCTION 

a  =  0.1  N  -  50  Cj  *  10.0  j  *  1, . . . ,50 
uj  *  (1.8  +  0.008j)x  ,j  s  1,...50 
At  «  0.05  n  =  1024 

A  typical  forcing  function  history  generated  by  using  these 
parameters  is  shown  in  Figure  8-12.  Actual  forcing  functions 
measured  in  the  field  usually  contain  a  certain  amount  of  noise. 
Inputs  with  noise  to  signal  ratios  of  six  and  eight  percent  are 
shown  in  Figures  8-13  and  8-14,  respectively. 

The  response  of  some  SDF  systems  to  the  forcing  input  were 
computed.  The  energy  dissipated  in  each  structure  is  listed  with 
the  structural  parameters  in  Tables  8.3A  and  8.38.  All  cases 
were  described  above.  The  notation  of  the  system  parameters  is 
as  follows:  k  Is  the  Initial  stiffness;  c  is  viscous  damping;  ky 
is  the  yield  stiffness;  zy  Is  the  yield  displacement;  z,^  is  the 
maximum  displacement  of  an  SDF  system;  m  is  the  mass  of  the  SDF 
structure. 


TABLE  8.3A.  SYSTEM  PARAMETERS 


TABLE  8.3B.  ENERGY  DISSIPATION  FOR  CASE  1  THROUGH  CASE  4 


Cases 

zy 

Damping  Energy 
Dissipated 

Spring  Energy 
Dissipated 

Total  Energy 
Dissipated 

Case  1 

00 

11028.20 

0.0 

11028.2 

Case  2 

6 

10199.30 

405.21 

10604.0 

Case  3 

5 

8364.27 

1186.74 

3551.0 

Case  4 

4 

6300.53 

1924.48 

8225.0 

First,  the  response  of  a  linear  system  (case  1)  was  com¬ 
puted.  Then,  a  slightly  nonlinear  response  of  an  SDF  structure 
was  computed  for  analysis  in  case  2.  The  displacement  response 
versus  time  for  case  2  is  plotted  in  Figure  8-15.  The  spring 
restoring  force  versus  displacement  is  shown  in  Figure  8-16.  A 
small  plastic  deformation  is  shown.  The  total  restoring  force 
versus  displacement  is  plotted  in  Figure  8-17.  The  measured 
responses  for  case  2  with  a  certain  amount  of  noise  are  plotted 
in  Figures  8-18  and  8-19.  Figure  8-18  shows  a  measured  signal 
with  six  percent  noise  to  signal  ratio.  Figure  8-19  shows  a 
measured  signal  with  ten  percent  noise  to  signal  ratio.  Two  more 
severe  nonlinear  responses  were  computed  for  analysis  in  cases  3 
and  4. 

The  displacement  response  versus  time  for  case  4  is  plotted 
in  Figure  8-20,  and  the  measured  response  for  case  4  with  ten 
percent  noise  to  signal  ratio  is  shown  in  Figure  8-21.  The 
spring  restoring  force  versus  displacement  is  shown  in  Figure 
8-22.  A  considerable  permanent  set  is  evident  in  Figure  8-20. 
Figure  8-22  shows  that  plastic  deformation  occurs  In  the  struc¬ 
ture  in  both  directions  of  motion. 


Restoring  Force  -  Kips  Restoring  Force  -  Kips 


[till 
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Using  the  forcing  function  input,  described  above,  and  the 
computed  responses,  the  parameters  of  the  structures  were  identi¬ 
fied.  The  results  of  the  parameter  identification  are  given  in 
Tables  8.4  through  8.7.  These  results  provide  the  identified 
parameters  in  the  noise-free  case.  The  energy  dissipated  by  the 


identified 

system 

is  listed  next 

to  the  identified 

parameters, 

For  method 

6,  the  parameters  ao. 

ai,  a2,  and  a3  Identify  with 

parameters 

k 0)  6, 

co,  and  a,  respectively. 

TABLE 

8.4. 

IDENTIFIED  PARAMETERS  AND 

ENERGY 

DISSIPATED 

FOR  CASE  1 

Method 

ao 

ai 

82 

83 

Energy 

1 

39.17 

1.26 

11210.0 

2 

39.17 

1.28 

0.0 

11090.0 

3 

40.01 

1.28 

10720.0 

4 

40.01 

2.48 

0.0236 

9270.0 

5 

39.50 

1.26 

0.0 

11200.0 

6 

39.33 

0.0 

1.274 

0.0 

10741.0 

TABLE 

8.5.  : 

IDENTIFIED  PARAMETERS  AND 

ENERGY 

DISSIPATED 

FOR  CASE  2 

Method 

ao 

ai 

a2 

a3 

Energy 

1 

33.27 

1.45 

7968.0 

2 

33.27 

3.98 

0.05 

7028.0 

3 

35.84 

1.28 

10500.0 

4 

35.84 

3.57 

0.042 

7690.0 

5 

34.27 

3.68 

0.062 

10370.0 

6 

36.56 

0.001 

1.094 

0.022 

10077.0 

Method 

ao 

ai 

32 

a3 

Energy 

1 

33.27 

1.67 

7376.0 

2 

33.27 

2.28 

0.014 

7053.0 

3 

34.61 

1.75 

7703.0 

4 

34.61 

2.50 

0.018 

7486.0 

5 

30.30 

3.95 

0.066 

7598.0 

6 

40.26 

-0.0018 

1.708 

0.001 

8186.5 

Section  4.0 

demonstrated 

that  when  the  higher  order 

linear  model  Is  used  to  simulate  the  actual  system  behavior,  the 
parameter  ao  must  be  estimated  first  In  the  identification  pro¬ 
cedures,  methods  1  through  4.  The  estimation  of  this  parameter 
can  be  executed  either  by  simply  searching  for  a  minimum  In 
Q(w)  ♦  e(w),  or  by  using  a  curve-fit  to  Q(w)  ♦  e(«).  and  then 
finding  the  minimum  of  the  curve.  Figure  8-23  shows  a  realiza¬ 
tion  of  Q(<*>)  for  a  specific  case.  This  Is  the  ratio  of  the 
Fourier  transform  moduli  of  the  structure  Input  and  response.  A 
example  of  the  quantity  Q(«)  +  e(«)  is  shown  In  Figure  8-24.  It 


QU) 


Frequency (w)  -  Second/cycle 
Figure  8-23.  A  realization  of  Q(u). 


Frequency(w)  -  Second/cycle 
Figure  8-24.  A  realization  of  Q(u>)  +  e(w) 
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is  apparent  from  this  diagram  why  the  use  of  a  curve-fit  provides 
better  results. 

Figures  8-25  through  8-27  show  comparisons  between  the 
responses  of  the  identified  systems  computed  by  different  meth¬ 
ods,  and  the  actual  response  of  the  bilinear  hysteretic  system  in 
case  2.  The  responses  of  the  identified  systems  match  the 
response  of  the  actual  system  so  closely  that  it  is  difficult  to 
distinguish  the  two  responses  in  these  figures.  More  identified 
responses  for  case  4  are  shown  in  Figures  8-28  through  8-30.  The 
model  responses  do  not  match  the  actual  response  as  closely  when 
residual  deformation  exists  In  the  actual  structure  since  the 
models  cannot  accumulate  permanent  deformation.  However,  peak 
responses  in  the  models  match  the  actual  system  response  quite 
well. 

In  this  example,  noise  was  added  to  the  forcing  function  and 
response  signals;  then  the  system  parameters  were  identified. 

The  results  are  summarized  In  Tables  8.8  and  8.9. 

» 

TABLE  8.8.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 

FOR  CASE  1  WITH  TEN  PERCENT  NOISE  TO  SIGNAL  RATIO 


Method 

ao 

ai 

a2 

33 

Energy 

1 

39.17 

1.37 

11022.0 

2 

39.17 

1.45 

0.003 

10900.0 

3 

39.93 

1.42 

10300.0 

4 

39.33 

2.44 

0.022 

9497.0 

5 

34.30 

2.98 

0.051 

11710.0 

6 

38.44 

0.0 

1.286 

0.0 

10929.0 

Displacement  -  In 


TABLE  8.9.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 


FOR  CASE  4 

WITH  TEN 

PERCENT  NOISE 

TO  SIGNAL  RATIO 

Method 

ao 

ai 

»2 

a3 

Energy 

1 

33.27 

1.89 

7092.0 

2 

33.27 

1.42 

0.021 

9980.0 

3 

32.02 

1.95 

6476.0 

4 

32.02 

2.60 

0.021 

6803.0 

5 

34.46 

2.70 

0.043 

11400.0 

6 

32.83 

0.0026 

1.543 

0.011 

7016.0 

These  results  show  that  the  parameter  identification  procedure  is 
still  effective  when  noise  is  present. 

8.3  Example  3 

In  this  example,  a  parameter  identification  problem  is 
solved  using  the  frequency  domain  approach.  The  methods  used  to 
identify  the  parameters  were  described  in  section  8.2,  namely 
methods  I  through  6.  The  same  forcing  function  as  illustrated  in 
Example  3  is  used.  The  only  difference  in  this  example  is  that 
the  response  was  simulated  by  a  second-order  time  varying  param¬ 
eter  system.  Unlike  the  responses  simulated  in  Example  2,  this 
example  is  a  linear  system  with  time  dependent  stiffness  and 
damping.  The  parameters  of  the  system  and  total  energy  dissi¬ 
pated  are  listed  in  Table  8.10. 

TABLE  8.10.  SYSTEM  PARAMETERS  FOR  CASE  5 

ko  =  39.48  co  *  1.257 

B  *  -  0.01  a  ■  0.01 

Total  Energy  «  9799.0 

The  definitions  of  the  symbols  are  the  same  as  In  Equation  5-2. 
This  case  was  described  in  section  8.2  as  case  5. 


The  displacement  response  versus  time  for  case  5  is  shown  in 
Figure  8-31.  The  total  restoring  force  versus  displacement  for 
this  case  is  illustrated  in  Figure  8-32.  Note  that  the  major 
axis  of  the  loops  depicted  in  the  diagram  have  different  slopes. 
This  occurs  because  the  system  stiffness  diminishes  with  time. 

The  parameters  identified  using  methods  1  through  6  together  with 
the  total  energy  dissipated  in  the  corresponding  systems  are 
shown  in  Table  8-11. 


TABLE  8.11.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 


FOR  CASE  5 

Method 

ao 

ai 

a2 

a3 

Energy 

1 

34.70 

1.36 

9266.0 

2 

34.70 

3.25 

0.039 

7883.0 

3 

37.15 

1.29 

10960.0 

4 

37.15 

2.82 

0.031 

9074.0 

5 

34.31 

2.97 

0.05 

11160.0 

6 

37.56 

-  0.006 

1.23 

0.019 

9120.2 

It  is  shown  that  methods  1  through  6  can  also  be  used  to  identify 
model  parameters  when  noise  is  present.  The  results  obtained 
when  the  measured  signals  contain  noise  are  shown  in  Table  8.12 


TABLE  8.12.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 

FOR  CASE  5  WITH  SIX  PERCENT  NOISE  TO  SIGNAL  RATIO 


Method 

ao 

ai 

a2 

as 

Energy 

1 

34.70 

1.38 

9419.0 

2 

34.70 

2.75 

0.03 

8623.0 

3 

37.34 

1.31 

11180.0 

4 

37.34 

2.24 

0.021 

10280.0 

5 

34.32 

2.49 

0.038 

11040.0 

6 

38.30 

-  0.007 

1.254 

0.025 

9245.6 
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Figure  8-30.  Comparison  between  actual  response 
(thin  line)  and  model  (method  6) 
response  (thick  line) 
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Figure  8-31.  Displacement  response  of  second-order 
time  varying  parameter  system  to  force 
in  Figure  8-15(Case  5) 
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Figures  8-33,  8-35,  and  8-36  compare  the  identified  system 
responses  to  the  response  of  the  actual  system.  The  simulated 
and  actual  responses  match  quite  well  in  all  cases.  The  model 
including  time  parameters  provides  the  best  match.  Figure  8-34 
shows  the  total  restoring  force  versus  displacement.  A  change  in 
slope  of  the  major  axis  of  the  loops  is  observed.  This  behavior 
matches  the  real  system  behavior  shown  in  Figure  8-32. 

8.4  Example  4 

In  this  example,  a  two  degree-of -freedom  (2DF)  system  is 
considered  and  its  parameters  are  identified  by  using  the  proce¬ 
dure  outlined  in  Chapter  6.  The  system  is  a  shear-beam  lumped 
mass  structure  as  shown  in  Figure  6-1. 

The  parameters  of  the  input  excitation  are  listed  in  Table 
8.13.  The  notation  for  the  parameters  is  the  same  as  in  previous 
examples.  A  typical  forcing  function  generated  using  these 
parameters  is  shown  in  Figure  8-37.  In  some  of  the  following 
analyses  noise  is  contained  in  the  measured  data.  An  input  with 
noise  to  signal  ratio  of  eight  percent  is  shown  in  Figure  8-38. 

TABLE  8.13.  PARAMETERS  OF  THE  FORCING  FUNCTION 

a  -  0.1  N  *  50  Cj  s  10.0  j  =  1,...,50 
At  =  0.05  n  =  1024  <*>j  =  (1.6  +  0.04  j)*,  j=l,...,50 

The  response  signals  used  to  represent  the  measured  signals 
were  generated  by  computer  program  BLNMDF.  BLNMDF  generated  the 
displacement,  velocity,  and  acceleration  response  of  a  bilinear 
MDF  system.  First,  the  response  of  a  linear  system  (case  1)  was 
computed.  Then,  a  nonlinear  response  was  computed  for  analysis 
in  case  2. 

Hie  structural  parameters  used  in  cases  1  and  2  in  this 
numerical  example  are  given  in  Table  8.14.  The  numerical  sub¬ 
scripts  on  the  parameters  refer  to  the  story  and  mass  numbers 
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Figure  8-32.  Total  restoring  force  versus  displace¬ 
ment  for  time  varying  parameter  system 
'Case  5) 
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Figure  8-33.  The  comparison  between  measured  response 
(light)  and  the  identified  response 
(dark)  by  method  6  for  Case  5 
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Figure  8-37. 


TABLE  8.14.  SYSTEM  PARAMETERS  (cases  1  and  2) 


ki 

=  39.48 

ci  *  1,257 

m  * 

1.0 

k2 

=  39.48 

c  2  *  1.257 

m  * 

1.0 

At 

=  0.05 

n  =  1024 

*di 

a  1.2,  y(j2  =  1.0  (for  case  2) 

in  the  2DF  system.  Since  yielding  only  occurs  in  case  2  the 
yield  level  parameters  are  only  used  in  case  2. 

The  displacement  response  versus  time  for  case  1  is  plotted 
in  Figures  8-39  and  8-40.  Note,  the  displacement  is  relative 
displacement  with  respect  to  each  degree-of -freedom.  The  rela¬ 
tive  displacement  histories  for  case  2  are  shown  in  Figures  8-41 
and  8-42.  Figures  8-43  and  8-44  show  the  measured  signals  with 
ten  percent  noise  to  signal  ratio  corresponding  to  the  actual 
response  measurements  in  Figures  8-41  and  8-42. 

The  energy  dissipated  by  the  structure  during  the  structural 
responses  of  cases  1  and  2  are  presented  in  Table  8.15.  Energy 
quantities  dissipated  in  both  first  and  second  stories  are 
given.  Note  that  similar  amounts  of  energy  are  dissipated  in  the 
linear  and  nonlinear  structures. 

TABLE  8.15.  ENERGY  DISSIPATION 

First  Story  Second  Story  Total 

Case  1  (Linear)  527.80  53.00  580.80 

Case  2  (Nonlinear)  447.20  105.10  552.30 

The  program  IDFID  performs  a  frequency  domain  parameter 
identification  for  the  second  and  third  order  linear  model.  It 
accepts  both  an  input  signal  from  FORCE  and  response  signals  from 
BLNM0F.  The  white  noise  signals  are  added  to  the  corresponding 
input  data  whenever  noises  are  included  in  the  analysis.  It 
should  be  noted  that  when  noise  signals  are  included  in  the 
responses,  a  different  signal  is  added  for  each  degree-of - 
freedom. 
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Figure  8-39.  Relative  displacement  response 
between  base  and  mass  1  in  2  DF 
structure.  Case  1  (linear) 


Figure  8-40.  Relative  displacement  response 
between  mass  1  and  mass  2  in  2DF 
structure.  Case  1  (linear) 
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Figure  8-43.  Relative  displacement  response 
between  base  and  mass  1  in  2DF 
structure  plus  10  percent  noise, 
(corresponding  to  Figure  8-41) 
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Figure  8-44.  Relative  displacement  response 
between  mass  1  and  mass  2  in  2DF 
structure  plus  15  percent  noise, 
(corresponding  to  Figure  8-42) 


MDFID  performs  the  necessary  Fourier  transforms  and  other 
data  operations.  Then  the  parameter  identification  is  executed. 
Once  the  parameters  are  known,  the  energy  dissipated  in  each 
degree-of -freedom  is  calculated. 

The  parameters  identified  using  the  program  MDFID  are  listed 
in  Tables  8.16  through  8.19.  The  system  parameters  identified 
for  case  1  (linear  response)  where  no  noise  is  included  in  the 
measured  signals,  are  given  in  Table  8.16.  The  system  parameters 
identified  for  case  1,  where  noise  is  included  in  the  measured 
signals,  are  given  in  Table  8.17. 

TABLE  8.16.  IDENTIFIED  PAAMETERS  AND  ENERGY  DISSIPATED 

FOR  CASE  1  WITHOUT  NOISE 
Convergence  Criteria  *  4  percent 

order  ki  k2  ci  C2  ai  a2  Ei  E2  E-p 

second  41.79  42.10  1.29  0.90  -  -  570.9  42.17  613.0 

third  41.30  43.21  1.76  0.605  0.0093  0.0014  634.5  29.66  664.2 

TABLE  8.17.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 

FOR  CASE  1  WITH  NOISE 
Convergence  Criteria  =  4  percent 

order  ki  k2  ci  C2  ai  a2  Ei  E2  E-p 

second  44.37  38.4  1.30  0.90  —  —  563.1  56.22  619.3 

third  43.45  39.07  1.38  0.96  0.0005  0.0037  592.7  42.80  635.5 

The  system  parameters  identified  for  case  2  (nonlinear) 
where  no  noise  is  Included  in  the  measured  signals  are  given  in 
Table  8.18.  The  system  parameters  identified  for  case  2,  where 
noise  is  Included  in  the  measured  signals,  are  given  in  Table 
8.19.  Each  table  lists  a  convergence  criterion  used  in  obtaining 
the  parameter  estimates.  These  quantities  reflect  a  limit  on  the 
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change  in  parameter  value  required  in  the  computation  before  the 
identification  computation  is  terminated.  The  smallest  conver¬ 
gence  criteria  values  yield  the  most  accurate  results. 

TABLE  8.18.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 

FOR  NONLINEAR  CASE  (CASE  2)  WITHOUT  NOISE 
Convergence  Criterion  =  4  percent 

order  ki  k2  ci  C2  ai  a2  Ei  E2  Ey 

second  29.04  33.87  1.81  1.81  —  -  861.3  81.43  948.73 

third  29.07  32.30  2.11  1.89  0.0082  0.0048  676.10  140.60  816.70 

TABLE  8.19.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 

FOR  CASE  2  WITH  NOISE 

Convergence  Criterion  *  2  percent 
order  ki  k2  ci  C2  ai  a2  Ei  E2  Ej 

second  35.33  19.91  1.21  1.38  -  -  434.8  142.6  577.3 

third  35.67  18.671  1.20  1.31  0.0033  0.001  397.1  133.2  530.3 

Comparisons  among  Tables  8.14,  8.16  and  8.17  show  that  the 
parameters  obtained  using  the  program  MDFID  provide  good  esti¬ 
mates  of  the  actual  system  parameters.  The  structural  responses 
predicted  using  the  second  and  third  order  models  are  shown  in 
Figures  8-45  through  8-48  where  they  are  compared  with  the  actual 
responses.  Model  responses  and  energy  dissipated  were  computed 
using  computer  program  HOMDF.  Figure  8-45  shows  the  relative 
displacement  response  between  the  base  and  mass  1  for  the  second 
order  model  (parameters  from  Table  8.16)  and  the  actual  response 
(Figure  8-39).  The  responses  are  almost  identical.  Figure  8-46 
shows  the  corresponding  relative  displacement  responses  between 
masses  1  and  2.  Figure  8-47  shows  the  relative  displacement 
response  between  the  base  and  mass  1  for  the  third  order  model 
(parameters  from  Table  8.16)  and  the  actual  response  (Figure 
8-39).  The  responses  are  almost  identical.  Figure  8-46  shows 
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Figure  8.45.  Comparison  of  second  order  model  and  actual 
relative  displacement  response  between  base 
and  mass  1.  (Case  1) 
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the  corresponding  relative  displacement  responses  between  masses 
1  and  2.  Figure  8-47  shows  the  relative  displacement  response 
between  the  base  and  mass  1  for  the  third  order  model  (parameters 
from  Table  8.16)  and  the  actual  response  (Figure  8-39).  Again, 
responses  are  almost  identical.  Figure  8-48  shows  the  corres¬ 
ponding  relative  displacement  responses  between  masses  1  and  2. 

The  energy  dissipated  during  structural  response  was  com¬ 
puted  using  the  second  and  third  order  models  in  program  HDMDF, 
using  the  parameters  identified  in  the  presence  of  no  noise  and 
the  parameters  identified  in  the  presence  of  noise.  The  results 
are  given  in  Table  8.16  and  8.17.  In  both  cases  good  agreement 
with  the  actual  energy  dissipated  is  found. 

Comparison  among  Tables  8.14,  8.18  and  8.19  show  that  the 
parameters  obtained  using  the  program  MDFID  are  changed.  This 
change  reflects  the  system  nonlinearity.  The  structural 
responses  predicted  using  the  second  and  third  order  models  for 
case  2  are  shown  in  Figures  8-49  through  8-52  where  they  are 
compared  with  the  actual  nonlinear  responses.  Model  responses 
and  energy  dissipated  were  computed  using  computer  program 
HDMDF.  Figures  8-49  and  8-50  show  the  relative  displacement 
responses  for  the  second  order  model  (parameters  from  Table  8.19) 
and  the  actual  responses  (Figures  8-41  and  8-42).  The  model 
responses  do  not  match  the  actual  response  as  closely  when  resid¬ 
ual  deformation  exists  in  the  actual  structure  since  the  models 
cannot  accumulate  permanent  deformation.  However,  peak  responses 
in  the  model  match  the  actual  system  response  quite  well.  Fig¬ 
ures  8-51  and  8-52  show  the  relative  displacement  responses  for 
the  third  order  model  (parameters  from  Table  8.19)  and  the  actual 
responses  (Figures  8-41  and  8-42).  Again,  responses  are  not 
closely  predicted  by  the  model.  However,  peak  responses  of  the 
model  match  the  actual  response  very  well. 
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Comparison  of  second  order  model  and  actual 
relative  displacement  response  between  base 
and  mass  1 .  (Case  2) 
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CHAPTER  9 

SUMMARY  AND  CONCLUSIONS 

The  objective  of  this  study  was  to  develop  approximate 
linear  models  for  the  simulation  of  inelastic  system  response  and 
the  measurement  of  damage  accumulation  in  a  structure.  It  was 
assumed  that  energy  dissipated  is  related  to  the  accumulation  of 
damage.  The  model  parameters  were  identified;  then  the  energy 
dissipated  during  a  strong  motion  was  calculated.  The  displace¬ 
ment  response  and  the  energy  dissipated  in  each  model  were  com¬ 
pared  with  the  displacement  response  and  energy  dissipated  in  the 
actual  structure. 

Three  basic  models  were  considered  in  this  study.  These  are 
second  and  third  order  linear  models  with  constant  coefficients, 
and  a  second  order  linear  model  with  time-varying  parameters. 

The  parameters  of  the  models  were  estimated  using  two  basic 
approaches,  namely  time  domain  and  frequency  domain  approaches. 
The  time  domain  approach  was  first  introduced  to  estimate  the 
model  parameters  by  using  the  least  squares  method  through  which 
the  modeling  error  is  minimized  with  respect  to  the  measured 
data.  One  reason  for  not  using  the  time  domain  approach  through¬ 
out  this  investigation  is  that  the  parameters  identified  using 
the  time  domain  approach  are  inaccurate  when  noise  signals  are 
present. 

Then,  the  frequency  domain  approach  was  used  to  identify  the 
model  parameters.  In  this  approach,  both  analytical  and  search 
techniques  were  applied  to  find  the  system  parameters.  A  good 
set  of  initial  guesses  of  the  parameters  was  an  important  concern 
when  the  search  technique  was  used  to  execute  the  identification 
process. 


128 


Several  numerical  examples  were  solved,  and  some  of  them  are 
summarized.  Experience  obtained  in  solving  the  numerical  exam¬ 
ples  lead  to  the  following  conclusions. 

(1)  Linear  and  nonlinear  hysteretic  SDF  systems  can,  in 
some  respects,  be  accurately  modeled  using  second-and 
third-order  linear  differential  equations  with  constant 
coefficients,  and  a  second-order  linear  differential 
equation  with  time-varying  coefficients.  Specifically, 
the  models  provide  accurate  simulation  when  displace¬ 
ment  response  and  energy  dissipated  criteria  are  used. 

(2)  A  direct,  time  domain  approach  can  be  used  to  identify 
model  parameters  when  the  force  input  and  acceleration 
response  measurements  are  not  noisy. 

(3)  The  frequency  domain  approach  can  be  used  to  identify 
model  parameters  of  all  three  models  when  the  force  and 
response  measurements  are  noisy. 

(4)  The  second  order  model  with  time-varying  coefficients 
provides  the  best  simulation  of  system  response  and 
energy  dissipated  among  the  three  models  considered. 

(5)  The  parameters  of  both  single-degree-of-freedom  (SDF) 
and  multi-degree-of-freedom  (MDF)  systems  can  be  iden¬ 
tified.  The  numerical  examples  show  that  one  and  two 
degree-of -freedom  systems  can  be  identified.  For  the 
model  presented  in  Chapter  6,  it  is  anticipated  that 
some  difficulty  would  exist  in  identifying  the  param¬ 
eters  of  a  system  with  three  or  more  degrees  of  freedom 
due  to  the  number  of  parameters  involved. 

(6)  The  energy  dissipated  in  a  structural  system  is  related 
to  system  damage. 


While  the  procedure  developed  in  this  investigation  provides 
means  for  the  simulation  of  response  and  the  estimation  of  damage 
in  inelastic  structures,  some  improvements  can  be  made.  The 
systems  considered  in  this  study  are  one  and  two-degree-of- 
freedom;  future  investigations  should  include  many  degree-of- 
freedom  structures.  The  models  used  in  this  study  do  not  permit 
the  accumulation  of  plastic  deformation;  future  investigations 
should  consider  models  that  allow  plastic  deformation  to  accumu¬ 
late.  The  tests  that  are  summarized  in  Chapter  7  show  that 
material  damage  is  related  to  energy  dissipated;  future  experi¬ 
ments  should  be  performed,  and  a  mathematical  model  characteriz¬ 
ing  the  results  should  be  developed.  Finally,  analyses  should  be 
performed  to  establish  the  spacial  distribution  of  energy  dissi¬ 
pated  in  actual  structural  members. 
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****************************************** 

*  PROGRAM  "BILIN"  * 

******************************************* 


SN:  EXCITATION  NOISE 
RN:  RESPONSE  NOISE 
MAXIMUN  DISPLACEMENT  EAQUL  TO  6.7 
INDEX=1  MEANS  NOISE  PRESENTED 

DIMENSION  SF( 1024) , DD( 1024) , SN( 1024) ,  RN( 1024) 
DIMENSION  AA(1024) ,RESF( 1024), SKI (1024) ,SPI( 1024) 
CALL  OPSYS ( ' ALLOC ' , ' RH1 ' , 10 ) 

CALL  OPSYS ( ' ALLOC ' , ' FD1 ' , 15 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH2 ' , 25 ) 

INDEX=1 
N=1024 
SM=1 . 0 
SC=1. 256637 
SK=39 . 48 
YSTI=0 . 0 
YDIS=4 . 0 
DT=0 . 05 

READ( 10, 50)  ( SF ( I ) , SN( I ) , 1=1 , N) 

READ (25, 55 )  (RN( I ) , I=1,N) 

CALL  BILIN( SF, SC, SM, SK, YSTI , YDIS,DT,N,DD, ENED, SPI ) 
IF( INDEX. EQ.l)  GO  TO  45 
DO  12  1=1, N 
T=DT* ( 1-1 ) 

RESF( I )=SF ( I )-SM*AA( I ) 

WRITE( 15, 50)T,DD( I ) 

12  CONTINUE 
GO  TO  60 
45  DO  40  1=1, N 
T=DT* ( I - 1 ) 

SF( I )=SF( I )+SN( I ) 

DD ( I ) =DD ( I ) +RN ( I ) 

40  WRITE ( 15 , 50)T, DD( I ) 

50  FORMAT (2E 12. 4) 

55  FORMAT (E12. 4) 

60  WRITE (6,*) ENED, ENED1,ENED2 
STOP 
END 

SUBROUTINE  FOR  GENERATING  THE  BILINEAR 
HYSTERETIC  SYSTEM. 

C7  :  DAMPING 

SM7 :  MASS 

SK7 :  STIFFNESS 

A7  :  YIELDING  STIFFNESS 

V7  :  YIELDING  DISPLACEMENT 

D9  :  TIME  INCREMENT 


V  :  OUTPUT  DISPLACEMENT 

VI  :  OUTPUT  VELOCITY 

VO  :  OFFSET  DISPLACEMENT 

V2  :  OUTPUT  ACCELERATION 

ENED: ENERGY  DISSIPATED 

SPI: SPRING  RESTORING  FORCE 

ENED1: TOTAL  SPRING  ENERGY  DISSIPATED 

ENED2: TOTAL  DAPPING  ENERGY  DISSIPATED 

SUBROUTINE  BILIN(F,  C7,  SM7,  SK7,  A7,  V7,D9,N,  V,ENED,  SPI ) 
DIMENSION  V(N) , VO ( 1024 ) , VI ( 1024 ) , V2 (1024) ,  F(N) 
DIMENSION  SPI (1024) 

U7=SK7*V7 
SK9=1 . 0-A7/SK7 

INITIALIZE  VARIABLES 

V(1)=0. 

V0(1)=0. 

V1(1)=0. 

V2(1)=F(1)/SM7 

START  THE  RESPONSE  CYCLE 

Q1=6.0*SM7/D9**2 
Q2=3 . 0*C7/D9 
Q3=6.0*SM7/D9 
Q4=3.0*SM7 
Q5=3 . 0*C7 
Q6=0.5*D9*C7 
Q7=3 . 0/D9 
Q8=3 . 0 
Q9=0.5*D9 
SK8=SK7 
NM=N-1 
ENED=0 . 0 
ENED 1=0 . 0 
ENED2=0 . 0 
DO  1199  1=1, NM 
11=1+1 

U1=Q1+Q2+SK8 
U2=Q3 *V1( I )+Q4*V2( I ) 

U3=Q5*V1 ( I ) +Q6*V2 ( I ) 

V5=(F( Il)-F( I )+U2+U3 )/Ul 
V6=Q7*V5-Q8*V1 ( I ) -Q9*V2 ( I ) 

V(I1)=V(I) +V5 
VI ( II )=V1 ( I ) +V6 

COMPUTE  THE  STIFFNESS  AT  T+DT 


X0=SK7  * ( V ( 1 1 ) -VO ( I ) ) 
X1=A7* ( V( II ) -V7 ) +U7 
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X2=A7* (V( II )+V7 )-U7 
I F ( XO . GT . XI ) GO  TO  1150 
I F ( X0 . LT . X2 ) GO  TO  1160 
SK8=SK7 
V0(I1)=V0(I) 

GO  TO  1170 

1150  IF  (V1(I1) .GT.0.0)SK8=A7 
IF  (V1(I1) .LE.0.0)SK8=SK7 
V0(I1)=(V(I1) -V7 ) *SK9 
GO  TO  1170 

1160  IF  (V1(I1) .LT.0.0)SK8=A7 
IF  (V1(I1) ,GE.O.O)SK8=SK7 
V0( I1)=SK9*(V{ Il)+V7) 

1170  SPI ( I1)=SK7* ( V( II) -V0( II ) ) 

V2(I1)=(F(I1)-C7*V1(I1)-SPI(I1) )/SM7 
ENED=ENED+D9*0.5*(V1(I)*(F(I)-SM7*V2(I) )+Vl(Il)* 
+(F(I1)-SM7*V2(I1) ) ) 

ENED1=ENED1+D9*0.5*(V1(I)*SPI(I)+V1(I1)*SPI(I1) ) 
ENED2=ENED2 +D9  *  0 . 5* ( VI ( I ) *C7*V1 ( I ) +V1 ( II ) *C7*V1 (II)) 
1199  CONTINUE 
RETURN 
END 


****************************************** 

*  PROGRAM  "  TIMEVA"  * 

****************************************** 


SOLVE  FOR  THE  SECOND  O.D.E. 

(TIME  VARIED  PARAMETERS) 

SN: INPUT  NOISE 
YN: RESPONSE  NOISE 

ALPHA: TIME  VARIED  COEFFICIENT  FOR  K 
BATA: TIME  VARIED  COEFFICIENT  FOR  C 

READ ( 5 ,  * )  CO, AKO, ALPHA, BATA 
CALL  OPSYS ( ' ALLOC ' , ' RH1 ' , 15 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH3 ' , 65 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH2 ' , 13 ) 

CALL  OPSYS ( ' ALLOC ' , 1 ID3 ' , 67 ) 

DIMENSION  DY (2 ) , Y(2 ) , F( 42 ) , SN( 1024) , YN( 1024) 
DIMENSION  SF( 1024), TS( 1024), S( 1024), INDEX( 1024) 
DIMENSION  Y3 ( 1024 ),V( 1024) 

EXTERNAL  SPLINE 

COMMON  /ARRAY/  TS, SF, S, INDEX 

NTS=1024 

DT=0 . 05 

READ (15, 8)  ( SF ( I ) , SN ( I ) , 1=1 , NTS ) 

READ ( 13 , 33 ) ( YN( I ) , 1=1, NTS) 

DO  2  1=1, NTS 
TS( I )=DT*( 1-1) 

2  CONTINUE 
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CALL  SPCOEF(NTS, TS, SF, S, INDEX) 

SM=1 . 0 
N=2 

TT=NTS*DT 
T=0 . 0 
Y( 1 )=0. 

Y ( 2 ) =0 . 0 
ENER=0 . 0 
10  L=3 
M=0 

50  CALL  RUNGE ( T , DT ,N,Y,DY,F,L,M,J) 

IF(M- 1)75, 10,75 
75  GO  TO  (100,200,999) ,L 
100  FQ=-FF(T) 

DY ( 1 ) =Y ( 2 ) 

DY(2 )=- ( 1 . +ALPHA*T) *AK0*Y( 1 )  -  ( 1 . +BATA*T) *C0*Y(2 ) +FQ 
GO  TO  50 
200  I=T/DT+1.1 

Y3  ( I )  =  ( 1 . +ALPHA*T) *AK0*Y ( 1 ) + ( 1 . +BATA*T ) *C0*Y( 2 ) 

V ( I ) =Y ( 2 ) 

IF( ( I— 1) . EQ.O)  GO  TO  111 

ENER=ENER+DT*0.5*( (Y3(I)*V(I) ) + (Y3 ( 1-1 ) *V( 1-1 ) ) ) 

GO  TO  250 

111  ENER=DT*0.5*(Y3(1)*V(1) ) 

NOISE  CASE 

112  SF( I )=SF( I )+SN( I ) 

YKK=Y ( 1 ) +YN ( I ) 

250  WRITE ( 65 , 800 )Y(1),Y3(I) 

IF(T-TT)260, 999 , 999 
260  GO  TO  50 
800  FORMAT (2E 12. 4) 

8  FORMAT (2E 12. 4) 

33  FORMAT (E12. 4) 

999  WRITE(6, * )  ENER 
STOP 
END 

REAL  FUNCTION  FF 
FUNCTION  FF(T) 

DIMENSION  TT( 1024) , FT( 1024) , SS( 1024) , INDE( 1024) 

COMMON  /ARRAY/  TT, FT, SS, INDE 

N=1024 

QA=SPLINE(N, TT, FT, SS, INDE,T) 

FF=QA 

RETURN 

END 

SOLVING  N-TH  ORDER  ODE 
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SUBROUTINE  RUNGE(T,DT,N, Y,DY, F,  L,M,  J) 
DIMENSION  DY (2),Y(2),F(42) 

GO  TO  ( 100, 110, 300) , L 

100  GO  TO  (101,110)  ,  IG 

101  J=1 
L=2 

DO  106  K  =1,N 

K1=K+3*N 

K2=K1+N 

K3=N+K 

F(K1)=Y(K) 

F(K3)=F(K1) 

106  F ( K2 ) =DY ( K ) 

GO  TO  406 

110  DO  140  K=1,N 
K1=K 
K2=K+5*N 
K3=K2+N 
K4=K+N 

GO  TO  (111,112,113,114),  J 

111  F(K1)=DY(K)*DT 
Y(K)=F(K4)+.5*F(K1) 

GO  TO  140 

112  F(K2)=DY(K)*DT 
GO  TO  124 

113  F(K3)=DY(K)*DT 
GO  TO  134 

114  Y(K)=F(K4)+(F(K1)+2.*(F(K2)+F(K3) ) +DY(K) *DT)/6 
GO  TO  140 

124  Y(K)=.5*F(K2) 

Y(K)=Y(K)+F(K4) 

GO  TO  140 

134  Y(K)=F(K4)+F(K3) 

140  CONTINUE 

GO  TO  (170,180,170,180), J 
170  T=T+.5*DT 
180  J=J+1 

IF( J-4)404, 404, 299 

299  M=1 

GO  TO  406 

300  IG=1 

GO  TO  405 

404  IG=2 

405  L=1 

406  RETURN 
END 


*  PROGRAM  "  URAND”  * 

*************** 


oooo  non  o  o  o  o  o  o  o  o  n  o  o  o  o  o  o  oo 
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THIS  PROGRAM  IS  FOR  GENERATING  THE  RANDOM  SIGNAL 

CALL  OPSYS ( ' ALLOC 1 , 1 RH2 ' , 66 ) 

CALL  OPSYS ( ' ALLOC ' , ’ RH1 '  ,  65 ) 

CALL  OPSYS ( ' ALLOC ' , *  RH4 '  ,11) 

REAL *8  URAND 

REAL  S ( 1024 ), SS ( 1024 ), DS ( 1024 ),T( 1024) ,TPF( 1024) 
REAL  EPS( 1024) ,SF( 1024) ,OUN( 1024) ,RENOI( 1024) 
DIMENSION  PS( 1024) , IND( 1024) , PHI ( 1024) 

DOUBLE  PRECISION  FD(1024) 

EXTERNAL  URAND 

INITIALIZATION 

WHEN  INDEX=1  SHOWS  NO  DERIVATIVE  FOR  FORCE, 
OTHERWISE  WILL  GO  TO  DERIVATIVE  INCLUDING  NOISE 

RATIO  : NO I SE/S I GNAL  RATIO 
S(I)  : FORCE 

DS(I)  : DERIVATIVE  OF  THE  FORCE 

SF ( I )  : FORCE  +  NOISE 

FD( I )  : DERIVATIVE  OF  ( FORCE+NOISE ) 

ZETA  : DAMPING  OF  THE  SYSTEM 
ALPHA  -.DECAY  RATIO 
NT  :  TOTAL  NUMBER  OF  POINT 

FN  : NATURAL  FREQUENCY 

NF  : NUMBER  OF  POINT  IN  THE  FREQUENCY  BAND 

INDEX=0 

RAT I 0=0 .0036 

ZETA— 0 . 1 

DT=0 . 05 

NF=50 

NT=1024 

NN=NT+1 

NSEED=0 

ALPHA=0 . 1 

C=10.0 

PI=4 . 0D0*DATAN( 1 . 0D0 ) 

FN=2.*PI 

CALCULATE  THE  FREQUENCY  AND  TIME 

F1=FN-(2.*ZETA*FN) 

IF(Fl.LT.O.O)  F1=0 . 0 
F2=FN+ ( 2 . *ZETA*FN) 

FB=F2-F1 
DF=(F2-F1 )/NF 
F3=F1- (DF/2 . 0 ) 

SI: VARIANCE  OF  THE  FORCE  NOISE 
S2: VARIANCE  OF  THE  OUTPUT  NOISE 


Sl=FB*RATIO*C*C/(2 . *DF) 


ooo  o  noon  nnon  n n o 


S2=PI *RATIO*C*C/( 8 . *DF* ( FN**3 ) *ZETA) 

DO  8  K=1,NF 
TPF ( K ) =F3  + ( K*DF ) 

PHI ( K ) =2 . *P I *URAND ( NSEED ) 

8  CONTINUE 

DO  9  J=1 , NT 
T( J)=(J-1)*DT 

9  CONTINUE 

GENERATE  THE  EXCITATION  AND  ITS  FIRST  DERIVATIVES 


DO  11  J=1 , NT 
S(1)=0. 

SS(1)=0. 

DO  11  K=1 , NF 

S(J)=S(J)+l.*C*COS(TPF(K)*T( J)-PHI(K) ) 

SS(J)=SS( J)+1.*C*TPF(K)*SIN(TPF(K)*T( J)-PHI(K) ) 

11  CONTINUE 

DO  12  J=1,NT 

S( J)=S( J)*EXP(-ALPHA*T( J) ) 

DS( J)=-SS(J)*EXP(-ALPHA*T(J) )-S(J)*ALPHA 

12  CONTINUE 

GENERATE  THE  RANDOM  EXCITATION  INCLUDES  NOISE. 

OUN: FORCE  NOISE, EPS:  RESPONSE  NOISE. 

CALL  NOISE ( OUN, NN, Si, 1) 

CALL  NOISE ( RENO I , NN , S2 , 2 ) 

CALL  NOI SE( EPS, NN,S2, NSEED) 

IF ( INDEX. EQ . 1 )  GO  TO  40 
FD( 1 )=0 . 0 
DO  15  1=2,  NT 

FD( I )=DS( I ) + (OUN( 1+1 ) -OUN( 1-1) )/(2. 0*DT) 

15  CONTINUE 

35  DO  38  J=1 , NT 

38  WRITE( 65, 22 )  S(J),DS(J) 

GO  TO  18 
40  DO  42  1=1, NT 

WRITE (77 , 21 )RENOI ( I ) 

WRITE(66,21)EPS( I ) 

42  WRITE ( 65 , 13 ) S( I ) , OUN( I ) 

13  FORMAT (2E12. 4) 

17  FORMAT(9X, 'TIME' ,11X, 'FORCE' , IX, 'FIRST  DERIVATIVE',/) 

21  FORMAT (1E12. 4) 

22  FORMAT (2E12. 4) 

18  STOP 
END 

SUBROUTINE  NOI SE( DPS, N, S2, NSEED) 


THIS  SUBROUTINE  GENERATES  A  SEQUENCE  OF  WHITE  NOISE 


ooo  oooo  oooo  oooo  oooo 
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IT'S  LENGTH  IS  N,  MEAN  ZERO,  VARIANCE  S2 

REAL* 8  URAND 
REAL  X( 12 ) ,DPS(N) 

C=SQRT(S2) 

DPS( 1 )=0 . 0 
DO  2  1=2, N 
DPS( I )=0. 0 
DO  1  J=l, 12 
X ( J ) =URAND ( NSEED ) 

1  DPS ( I )=DPS ( I ) +X(  J ) 

2  DPS( I )=(DPS ( I ) -6 . 0) *C 
RETURN 
END 

REAL  FUNCTION  URAND*8(IY) 

INTEGER  IY 

URAND  IS  A  UNIFORM  NUMBER  GENERATOR  BASED  ON  THE 
THEORY  AND  SUGGESTIONS  GIVEN  IN  D.E.  KNUTH  (1969). 

INTEGER  IA, IC, ITWO, M2 , M, MIC 
DOUBLE  PRECISION  HALFM 
REAL  S 

DOUBLE  PREC  x  S I  ON  DATAN , DSQRT 
DATA  M2/0/, I TWO/2/ 

IF  (M2  .NE.  0)  GO  TO  20 

IF  FIRST  ENTRY,  COMPUTE  MACHINE  INTEGER 
WORD  LENGTH. 

M  =  1 
10  M2  =  M 

M  =  ITWO*M2 

IF  (M  .GT.  M2)  GO  TO  10 
HALFM  =  M2 

COMPUTE  MULTIPLIER  AND  INCREMENT  FOR 
LINEAR  CONGRUENT I AL  METHOD 

IA  =  8*IDINT(HALFM*DATAN(1.D0)/8.D0)  +  5 
IC  =  2  * IDINT ( HALFM* ( 0 . 5D0-DSQRT ( 3 . DO ) /6 . DO ) )  +  1 
MIC  =  (M2  -  IC)  +  M2 

S  IS  THE  SCALE  FACTOR  FOR  CONVERTING  TO 
FLOATING  POINT. 

S  =  0 . 5 /HALFM 

COMPUTE  NEXT  RANDOM  NUMBER 

20  IY  =  IY*IA 
C 


THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHICH 
DO  NOT  ALLOW  INTEGER  OVERFLOW  ON  ADDITION 


IF  (IY  .GT.  MIC)  IY  =  (IY  -  M2)  -  M2 
IY  =  IY  +  IC 

THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHERE 
THE  WORD  LENGTH  FOR  ADDITION  IS  GREATER  THAN 
FOR  MULTIPLICATION 

IF  ( IY/2  .GT.  M2)  IY  =  (IY  -  M2)  -  M2 

THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHERE 
INTEGER  OVERFLOW  AFFECTS  THE  SIGN  BIT 

IF  (IY  .LT.  0)  IY  =  (IY  +  M2)  +  M2 

URAND  =  FLOAT ( IY)  *S 

RETURN 

END 


****************************************** 

*  PROGRAM  "  BLNMDF"  * 

****************************************** 


THIS  PROGRAM  GENERATE  THE  MDF  BILINEAR 
HYSTERETIC  STRUCTURE  RESPONSE 


N  : 
NT  : 
DT  : 
MC( I )  : 
CC( I )  : 
KC(I)  : 
AC  ( I )  : 
ZOI(I)  : 
ZOF(I) : 
Z1I(I): 
Z1F(I) : 
Z2F ( I ) : 
PSF(I)  : 


NUMBER  OF  DEGREE  OF  FREEDOM 
NUMBER  OF  TIME  STEP 
TIME  INCREMENT 

MASS  FOR  EACH  DEGREE  OF  FREEDOM 

DAMPING  FOR  EACH  DEGREE  OF  FREEDOM 

STIFFNESS  FOR  THE  SYSTEMS 

YIELD  STIFFNESS 

INITIAL  DISPLACEMENT 

FINAL  DISPLACEMENT 

INITIAL  VELOCITY 

FINAL  VELOCITY 

FINAL  VELOCITY 

FINAL  PERMANENT  DISPLACEMENT 


CALL  OPSYS ( ' ALLOC ' , ' MDF1 ' , 10 ) 

CALL  OPSYS ( ' ALLOC ' , ' MDF2 ' , 20 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH1 ' , 23 ) 

CALL  OPSYS( 'ALLOC' , ' RH2 ' ,25) 

CALL  OPSYS ( ' ALLOC ' , ’ RH3 ' , 40 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH4 ' , 46 ) 

DIMENSION  DUMV1 ( 2 ) , DUMV2 ( 2 ) , DUMM1 (2,2), DUMM2 (2,2) 
DIMENSION  DLO ( 2 ) , DL1 ( 2 ) , RU ( 2 ) , RL (2), PS 1(2), YD ( 2 ) 
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DIMENSION  CC(2) ,AC(2),Q(2) ,PSF(2) ,FF( 1024), RN( 1024) 
DIMENSION  CM(2 , 2 ) , RV( 2 ) , DZ0(2 ) , DZ1 ( 2 )  , DZ2 ( 2 ) 
DIMENSION  ZO 1(2), 211(2), Z2 1(2), F( 1024,2 ) , Z0F{2 ) 
DIMENSION  A1 ( 2 , 2 ) , A2 ( 2 , 2 ) , A3 ( 2 , 2 ) , A4 ( 2 ) , A5 ( 2 ) 
DIMENSION  AKMC ( 2 ) , AMC ( 2 ) , AKC ( 2 ) , AMM (2,2), AKM (2,2) 
DIMENSION  A1 1 ( 2 , 2 ) , A2 2 ( 2 ) , ZM ( 2 ) , ZS ( 2 ) , ZN ( 2 ) , ZT ( 2 ) 
DIMENSION  All I (2,2), ENED1 ( 1024) , ENED2 ( 1024) , DN( 1024) 
DIMENSION  Z1F(2 ) , Z2F(2 ) , FN( 1024) 

C 

C  INOISE=l  MEANS  NOISE  INCLUDED 
C 

INOISE=0 

N=2 

NT=1024 
DT=0 . 05 
AMC( 1 )=1 . 0 
AMC (2 )=1 . 0 
CC( 1 )=1 . 256 
CC(2 )=1 . 256 
AKC( 1 )=39 . 48 
AKC(2 )=39 . 48 
AC ( 1 ) =0 . 0 
AC ( 2 ) =0 . 0 
Z0I(1)=0.0 
ZOI (2 )=0 . 0 
Z1I(1)=0.0 
Z1I(.2)=0.0 
YD ( 1 )=1 . 2 
YD(2 )=1 . 0 
Nl=NT-l 
ENED1 ( 1 ) =0 . 0 
ENED2 ( 1 ) =0 . 0 
C 

C  READ  THE  FORCE  VECTOR 

C  FN(I): INPUT  NOISE 

C  RN(I) : RESPONSE  NOISE 

C 

READ(23 , 33 ) ( FF ( I ) , FN( I ) , 1=1 , NT) 

READ (25, 35 ) (RN( I ) , 1=1, NT) 

READ (46, 35 ) (DN( I ) , 1=1 , NT) 

33  FORMAT (2E12. 4) 

35  FORMAT (E12. 4) 

DO  30  1=1, NT 
F( I , 1 )=0 . 0 
30  F( I , 2 )=FF( I ) 

C 

C  FORM  THE  MM  AND  CM  MATRICES 
C 

DO  1111  1=1, N 
DO  1111  J=1 , N 
AKM( I , J)=0. 

AMM( I , J)=0 . 


oo  ooo  no o  noon  non  non 
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1111  CM ( I , J ) =0 . 

FORM  MASS  MATRIX 


DO  1112  1=1, N 
1112  AMM (1,1) =AMC ( I ) 

FORM  DAMPING  MATRIX 


DO  1113  1=1, N 
IP=I+1 

CM (1,1) =CC ( I ) 

I F ( I . EQ . N ) GO  TO  1113 
CM( I , I )=CM( 1,1) +CC( IP ) 

CM( I , IP )=-CC( IP ) 

CM( IP, I )=-CC( IP) 

1113  CONTINUE 

FORM  CONSTANTS  VECTORS  AND  MATRICES 
FOR  ITERATION. 

CDl=DT/2 . 0 
CD2=6 . 0/DT 
CD3=6 . 0/ ( DT*  *2 ) 

CD4=3 . 0/DT 

Dl=DT/2 . 0 

D2=DT*DT/6 . 0 

D3=DT*DT/2 . 0 

DO  1119  1=1, N 

A4 ( I ) =AC ( I ) / ( AKC ( I ) -AC ( I ) ) 

A5 ( I ) = ( AKC ( I ) -AC ( I ) ) /AKC ( I ) 

DO  1119  J=1,N 

A1 ( I , J)=AMM( I , J) +D1*CM( I , J) 

A2  ( I , J ) =DT*CM ( I , J ) 

1119  CONTINUE 


INVERSE  THE  MASS  MATRIX 
CALL  MATINV( AMM, A3 ,N) 

INITIATE  KM  AND  RESPONSE  VECTORS 


DO  1130  1=1, N 
IP=I+1 

AKM( I , I )=AKC( I ) 

IF( I . EQ . 2 )  GO  TO  1130 
AKM( I , I )=AKM( 1,1) +AKC( IP) 
AKM( I , IP )=-AKC( IP) 

AKM( IP, I )=-AKC( IP) 

1130  CONTINUE 


SET  INITIAL  CONDITIONS 


DO  1132  1=1 , N 
Z0I(I)=0. 

Z1I(I)=0. 

1132  PSI(I)=0. 

DO  1133  1=1, N 

1133  DUMV1 ( I ) =F (1,1) 

FIND  INITIAL  ACCELERATION  Z2I 

CALL  MVM ( A3 , DUMV1 , Z2 I , N ) 

WRITE  INITIAL  VALUES 

WRITE (6, 3) 

T=0 . 0 

WRITE(40, 100)  FF( 1 ) 

WRITE (10, 14) (ZOI(I) , 1=1 ,N) 

WRITE (20, 14) (PSI(I), 1=1, N) 

START  THE  RESPONSE  COMPUTATION 

DO  8999  IND=1,N1 
T=IND*DT 

COMPUTE  DZ2 

DO  1121  I=1,N 
DO  1121  J=1 , N 

1121  DUMM1 ( I , J ) =A1 ( I , J ) +D2  * AKM ( I , J ) 

CALL  MATINV(DUMM1,DUMM2,N) 

DO  1212  1=1, N 
DO  1212  J=1 , N 

1212  DUMM1 ( I , J ) =A2 ( I , J ) +D3  * AKM ( I , J ) 

CALL  MVM ( DUMM1 , Z2 I , DUMV1 , N ) 

DO  1213  1=1, N 
DO  1213  J=1 , N 

1213  DUMM1 ( I , J ) =DT* AKM ( I , J ) 

CALL  MVM(DUMM1 , Z1I , DUMV2 , N) 

DO  1214  1=1, N 

1214  DUMV1 ( I ) =- ( DUMV1 ( I ) +DUMV2 ( I ) ) + 

+(F( IND+1, I)-F(IND,I)) 

CALL  MVM ( DUMM2 , DUMV1 , DZ2 , N ) 

COMPUTE  DZ1  AND  DZO 

DO  1221  1=1 , N 
DZ1 ( I )=DT*Z2I ( I ) +DZ2 ( I ) *D1 
122 1  DZO ( I ) =D3*Z2 1(1) +D2  *DZ2 ( I ) +DT*Z1 1 ( I ) 


COMPUTE  Z1F  AND  ZOF 


noo  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
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DO  100  1=1 ,N 
DO  100  J=1,N 

100  All ( I , J )=AKM ( I ,  J ) +CD3 *AMM( I , J) +CD4*CM( I ,  J ) 
DO  110  I=1,N 

ZM( I )=CD2*Z1I ( I ) +3 . 0*Z2I ( I ) 

ZS(I)=3.0*Z1I(I )+CDl*Z2I ( I ) 

110  CONTINUE 

CALL  MVM ( AMM ,  ZM ,  ZN ,  N ) 

CALL  MVM ( CM ,  ZS , ZT , N ) 

DO  120  1=1 , N 

120  A22 ( I )=F ( IND+1 , I )  —  F ( IND, I ) +ZN( I ) +ZT( I ) 

CALL  MATINV(A11, A11I ,N) 

CALL  MVM (All I ,  A22 , DZO , N) 

DO  140  1=1  ,N 

140  DZ1 ( I )=CD4*DZ0 (I)-3.0*Z1I(I ) -CD1*Z2 1(1) 

DO  1231  1=1 ,  N 

Z1F ( I ) =Z1 1(1) +DZ1 ( I ) 

1231  ZOF ( I )=Z0I ( I ) +DZ0 ( I ) 

COMPUTE  THE  END  OF  STEP  KM  MATRIX 

DLO ( 1 )=Z0F ( 1 ) 

DL1 ( 1 )=Z1F ( 1 ) 

I F ( N . EQ . 1 ) GO  TO  1245 

DO  1241  I=2,N 

DLO ( I )=Z0F ( I )-Z0F( 1-1) 

1241  DL1 ( I )=Z1F( I )-ZlF( 1-1 ) 

1245  CONTINUE 

DO  1251  1=1, N 

RU ( I )=AKC( I ) * ( A4( I ) *PSI ( I ) +YD( I ) ) 

1251  RL( I )=AKC( I ) * ( A4( I ) *PSI ( I ) -YD( I ) ) 

DO  1272  1=1 ,N 

QCC=AKC( I ) * (DLO ( I ) -PSI ( I ) ) 

I F ( QCC . GT . RU ( I ) ) GO  TO  1258 
I F ( QCC . LT . RL ( I ) ) GO  TO  1265 
AKMC( I )=AKC( I ) 

PSF ( I )=PSI ( I ) 

GO  TO  1271 

1258  I F ( DL1 ( I ) . LE . 0 . ) GO  TO  1260 
AKMC ( I ) =AC ( I ) 

GO  TO  1261 

1260  AKMC ( I ) =AKC ( I ) 

1261  PSF ( I ) =A5 ( I ) * (DLO ( I ) -YD ( I ) ) 

GO  TO  1271 

1265  IF(DL1 ( I ) . GE . 0 . )GO  TO  1267 
AKMC ( I ) =AC ( I ) 

GO  TO  1268 

1267  AKMC ( I ) =AKC ( I ) 

1268  PSF ( I ) =A5 ( I ) * ( DLO ( I ) + YD ( I ) ) 

1271  Q( I ) =AKC ( I ) * (DLO ( I ) -PSF ( I ) ) 

1272  CONTINUE 

DO  1281  J=1 , N 


ooo  nooo  ooo  ooo 


JP=J+1 

AKM ( J ,  J ) =AKMC ( J ) 

IF( J.EQ.N)GO  TO  1281 

AKM  ( J ,  J )  =AKM  ( J ,  J )  + AKMC  ( JP ) 

AKM ( J , JP ) =- AKMC ( JP ) 

AKM ( JP , J )=- AKMC ( JP ) 

1281  CONTINUE 


COMPUTE  THE  RESTORING  FORCE  VECTOR 

DO  1291  1=1,  N 
RV ( I ) =Q ( I ) 

IF( I . LT.N)RV( I )=RV( I)-Q(I+1) 

1291  CONTINUE 

COMPUTE  Z2F 

CALL  MVM(CM, Z1F, DUMV1, N) 

DO  1311  1=1, N 

1311  DUMV2 ( I )=F( IND+1 , I ) -DUMV1 ( I ) -RV ( I ) 

CALL  MVM(A3,DUMV2,Z2F,N) 

PRINT  THE  RESULTS 

WRITE (6,3) 

IF( INOISE.NE. 1)  GO  TO  555 
ZNO I SE=ZOF ( 1 ) +RN ( IND+ 1 ) 

YNOISE=ZOF ( 2 ) +DN( IND+1 ) 

WRITE ( 10, 14)ZN0ISE, YNOISE 
YF=FF ( IND+ 1 ) +FN ( IND+ 1 ) 

WRITE ( 40, 100) YF 
GO  TO  558 

555  WRITE (10,14) ZOF ( 1 ) , ZOF ( 2 ) 

WRITE ( 40 , 100 ) F ( IND+ 1,2) 

558  WRITE (20, 14) (PSF( I ) , I=1,N) 

IP=IND+1 

DIS=Z0F(2 ) -Z0F( 1 ) 

REV=Z1F(2 ) -Z1F( 1 ) 

REACC=Z2F(2 ) -Z2F( 1 ) 

ENED1 ( IP)=Z1F (1)*(CC(1)*Z1F(1) +AKC( 1 ) *Z0F ( 1 ) ) 
ENED2 ( IP )=REV* ( CC ( 2 ) +REV+AKC ( 2 ) *DI S ) 

INITIATE  THE  NEXT  COMPUTATION  CYCLE 

DO  1341  1=1, N 
ZOI ( I )=Z0F( I ) 

Z1I ( I )=Z1F( I ) 

Z2I ( I )=Z2F( I ) 

PSI ( I )=PSF ( I ) 

1341  CONTINUE 
8999  CONTINUE 

CALL  SIMP(ENED1 , ENER1 , DT, NT) 


non  o  o  o 


CALL  S IMP ( ENED2 , ENER2 , DT , NT ) 
ENER=ENER1 +ENER2 
WRITE (6, 17)  ENER1 , ENER2 , ENER 
3  FORMAT ( / ) 

14  FORMAT (2E12. 4) 

17  FORMAT (3E12. 4) 

100  FORMAT (1E12. 4) 

STOP 

END 

FIND  THE  AREA  UNDER  THE  CURVE 

SUBROUTINE  SIMP(X,E,D,N) 
DIMENSION  X( 1024) 

N2=N/2-l 
E=X ( 1 ) 

DO  1  1=1, N2 
IA=2*I 
IB=IA+1 
E=E+4.0*X(IA) 

E=E+2 . 0*X( IB) 

1  CONTINUE 
N3=N-1 

E=E+0 . 5*X(N3 ) +1 . 5*X(N) 

E=D*E/3 . 0 

RETURN 

END 

SUBROUTINE  MULTIPLICATION 

SUBROUTINE  MVM(A,B,C,N) 
DIMENSION  A(2/2),B(2),C(2) 

DO  1  1=1 ,  N 
C(I)=0. 

DO  1  J=1 , N 

C( I )=A( I ,  J ) *B( J) +C( I ) 

1  CONTINUE 
RETURN 
END 


1 


i 


i 
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****************************************** 

*  PROGRAM  "TIMID''  * 

****************************************** 

MAIN  PROGRAM  FOR  LEAST- SQUARE  STRUCTURE 
IDENTIFICATION 

EXTERNAL  SUBS 
EXTERNAL  MULTI 

DIMENSION  C(720 , 2 ) , D( 10, 10) , AT(2, 720) ,CM(2,2) 
DIMENSION  F( 720, 1 ) ,  AA(2 ,  1 ) , CA( 720, 1 ) , ERR ( 720, 1 ) 
DIMENSION  CMS (2, 2) ,DIA(2,1) , CMS1 (2 , 2 ) , UM(2 , 2 ) 
DIMENSION  CM1 (2, 2), CT (2, 720), ERRT( 1,720) , AJN(2 , 2 ) 
DOUBLE  PRECISION  AT, CM, CT, AA, C, F, CA, ERR, ERRT 
DOUBLE  PRECISION  AJN, CM1, CMS,DIA, CMS1,UM 
CALL  OPSYS( 'ALLOC' ID1' ,33) 

DATA  FILE  NEED  TO  CHANGE  ACCORDING  THE 
SYSTEM  USED 

N  NUMBER  OF  POINTS  FOR  THE  IDETI FI CATION 
IM  NUMBER  OF  COLUM 

CALL  OPSYS ( ' ALLOC ' , ' ID1 ' , 34 ) 

READ ( 5 , * )  N,M, IM 
DO  10  1=1, N 

10  READ( 33 , 15 )  (C(I,J),  J=1,IM) 

DO  11  1=1,  N 

11  READ( 34, 16 )  (F(I,1)) 

CALL  TRAN(C,N, IM, AT) 

CALL  MULTI (AT, IM, N, C, IM, CMS1 ) 

DO  8  1=1, IM 
DO  8  J=1 , IM 
8  CM ( I , J ) =CMS1 ( I , J ) 

WRITE( 6, 30 ) 

DO  3  1=1, IM 

3  WRITE( 6, 4)  (CMS1 (I,J),J=1, IM) 

CALL  INVER(CMS1, IM,D,M, DETER) 

CALL  MULTI (CM, IM, IM,CMS1, IM,UM) 

WRITE ( 6, 33 ) 

DO  5  1=1, IM 

5  WRITE{ 6, 7 )  (UM(I, J) , J=l, IM) 

CALL  MULTI (CMS1, IM, IM, AT, N,CT) 

CALL  MULTI (CT, IM, N, F, 1 , AA) 

WRITE( 6, 34) 

DO  22  1=1, IM 
22  WRITE( 6,12)  AA(I,1) 

CALL  MULTI (C, N, IM, AA, 1 , CA) 

CALL  SUBS (F,N,1,CA, ERR) 

CALL  TRAN ( ERR, N, 1,ERRT) 

CALL  MULTI (ERRT, 1, N, ERR, 1, VAR) 

VAR=VAR/(N-IM) 

DEV=SQRT ( VAR ) 


WRITE( 6, 17 ) 

WRITE( 6, 37 )  DEV 
4  FORMAT (2E16. 6) 

7  FORMAT (2E16. 6) 

12  FORMAT (E16. 6) 

15  FORMAT (2E 12. 4) 

16  FORMAT ( 1E16 . 6 ) 

17  FORMAT  (5X, 'THE  DEVIATION  IS',/) 

20  FORMAT (6E12. 4) 

24  FORMAT (6E 12. 4) 

30  FORMAT (5X, 'MATRIX  Z',/) 

31  FORMAT(5X/ ' INVERSE  MATRIX',/) 

32  FORMAT(5X, ' INVERSE  MATRIX  IMPROVED',/) 

33  FORMAT (5X, 'UNIT  MATRIX  FOR  CHECKING',/) 

34  FORMAT ( 5X, ' INDENT I FIED  PARAMETERS',/) 

35  FORMAT(5X, 'MATRIX  CT',/) 

37  FORMAT (E18. 6) 

STOP 

END 

THIS  SUBROUTINE  IS  FOR  MATRIX  INVERSION 
SUBROUTINE  I NVER ( A , N , B , M , DET ) 

DOUBLE  PRECISION  A(2 , 2 ) , B(2 , 2 ) , IPVOT( 2 ) , INDEX(2 , 10 ) 
DOUBLE  PRECISION  T,PIVOT(2) 

COMMON  IPVOT, INDEX, PIVOT 
EQUIVALENCE  ( I ROW , JROW ) , ( I COL , JCOL ) 

INITIALIZATION 

DET=1 . 0D0 
DO  7  J=1 , N 
7  IPVOT (J)=0 
DO  135  1=1, N 

SEARCH  FOR  PIVOT  ELEMENT 


T=0 . ODO 


DO  9  J=1 , N 

I F ( I P VOT ( J ) - 1 )  13,9,13 
13  DO  23  K=1,N 

IF( IPVOT(K) -1 )43 ,23,81 
43  IF(DABS(T)-DABS(A( J,K) ) )  83,23,23 
83  IROW=J 


ICOL=K 


T=A ( J , K ) 

23  CONTINUE 
9  CONTINUE 

I P VOT ( I COL ) = I P VOT ( I COL ) + 1 


PUT  PIVOT  ELEMENT  ON  DIAGONAL 


I F ( IROW- ICOL)  73,  109,73 
73  DET=-DET 
DO  12  L=1 , N 
T=A( IROW, L) 

A( IROW, L)=A( ICOL, L) 

12  A( ICOL, L)=T 

IF(M)  109,109,33 
33  DO  2  L=1,M 
T=B ( I ROW , L ) 

B( IROW, L)=B( ICOL,  L) 

2  B( ICOL, L)=T 
109  INDEX ( I , 1 )=IROW 
INDEX( I , 2 )=ICOL 
PIVOT( I )=A( ICOL, ICOL) 
DET=DET*PIVOT( I ) 

DIVIT  PIVOT  ROW  BY  PIVOT  ELEMENT 

A( ICOL, ICOL)=l. 

DO  205  L=1,N 

205  A( ICOL, L)=A( ICOL, L)/PIVOT( I ) 
IF(M)  347,347,66 
66  DO  52  L=1,M 

52  B( ICOL, L)=B( ICOL, L)/PIVOT( I ) 

REDUCE  NON-PIVOT  ROW 

347  DO  135  LI=1,N 

IF(LI-ICOL)  21,135,21 
21  T=A(LI , ICOL) 

A(LI, ICOL)=O.ODO 
DO  89  L=1,N 

89  A (LI , L)=A(LI , L) -A( ICOL, L) *T 
IF(M)  135,135,18 
18  DO  68  L=1,M 

68  B(LI , L)=B( LI , L) -B( ICOL, L) *T 
135  CONTINUE 

INTERCHANGE  COLUMNS 


222  DO  3  1=1, N 
L=N- I + 1 

IF( INDEX(L, 1 ) - INDEX ( L, 2 ) )  19,3,19 
19  JROW= INDEX ( L , 1 ) 

JCOL=INDEX( L, 2 ) 

DO  549  K=1,N 
T=A(K, JROW) 

A(K, JROW)=A(K, JCOL) 

A(K, JCOL)=T 
549  CONTINUE 
3  CONTINUE 
81  RETURN 


THIS  SUBROUTINE  IS  FOR  MATRIX  MULTIPLICATION 


SUBROUTINE  MULTI ( A, N, M, B , L, C ) 

DIMENSION  A(N,M),B(M,L),C(N,L) 

DOUBLE  PRECISION  A,B,C 
DO  5  1=1, N 
DO  5  J=1,L 
C(I,  J)=O.ODO 
DO  5  K=1,M 

5  C(I, J)=C(I, J)+A(I,K)*B(K, J) 

RETURN 

END 

SUBROUTINE  MULTIA(TA, N,M, TB, L, TC) 

DIMENSION  TA(N, M) , TB (M, L) , TC(N, L) 

DOUBLE  PRECISION  TA , TB , TC 
DO  15  1=1, N 
DO  15  J=l, L 
TC ( I , J ) =0 . 0D0 
DO  15  K=1 , M 

15  TC( I , J)=TC( I , J) +TA( I , K) *TB(K, J) 

RETURN 

END 

THIS  SUBROUTINE  IS  FOR  MATRIX  TRANSFORMATION 

SUBROUTINE  TRAN( A, N, IM, AT) 

DIMENSION  A(N, IM) ,AT(IM,N) 

DOUBLE  PRECISION  A, AT 
DO  3  1=1, N 
DO  3  J=1 , IM 
3  AT( J, I )=A( I , J) 

RETURN 

END 

THIS  SUBROUTINE  IS  FOR  MATRIX  SUBSTRACTION 

SUBROUTINE  SUBS( A,N, IN, B, C) 

DOUBLE  PRECISION  A(N, IN) , B(N, IN) , C(N, IN) 

DO  101  1=1, N 
DO  101  J=1 , IN 
101  C( I , J )=A( I,J)-B(I,J) 

RETURN 

END 

SUBROUTINE  SET( A, IM, B, DIA) 


DIMENSION  A( IM, IM),B(IM, IM) , DIA( IM, 1 ) 
DOUBLE  PRECISION  A, B,DIA 
DO  3  J=l, IM 
DO  3  1=1, IM 


B(I,J)=A(I,J)/A(J,J) 

DIA( J, 1)=A( J, J) 

3  CONTINUE 
RETURN 
END 

THIS  SUBROUTINE  IS  FOR  IMPROVING  THE  ACCURACY 
OF  THE  MATRIX  INVERSION. 

SUBROUTINE  RIVISE(A,A1, IM,AIN) 

DOUBLE  PRECISION  A{2 , 2 ) , A1 (2 , 2 ) , AIN(2 , 2 ) # AJN(2 , 2 ) 
DOUBLE  PRECISION  UN ( 2 , 2  ) , UNI (2,2), SAI (2,2) 

DOUBLE  PRECISION  ERR (2,2), AA1 (2,2) 

UNI ( 1, 1)=1. 0D0 
UN1(2,2)=1.0D0 
UNI ( 3 , 3 ) =1 . 0D0 
UNI (4, 4)=1 . 0D0 
UNI ( 1 , 2 ) =0 . 0D0 
UNI ( 1 , 3 ) =0 . 0D0 
UNI ( 1 , 4 ) =0 . ODO 
UNI (2,1) =0 . ODO 
UNI ( 2 , 3 ) =0 . ODO 
UNI ( 2 , 4 ) =0 . ODO 
UNI ( 3 , 1 ) =0 . ODO 
UNI ( 3 , 2 ) =0 . ODO 
UNI (3,4) =0 . ODO 
UNI ( 4 , 1 ) =0 . ODO 
UNI ( 4 , 3 ) =0 . ODO 
DO  30  1=1, IM 
DO  30  J=l,  IM 

30  UN ( I , J ) =UN1 ( I , J ) *2 . ODO 
IC0UNT=1 

2  CALL  MULTI (A, IM, IM, A1 , IM, AA1 ) 

CALL  SUBS (UN, IM, IM, AA1 , SAI ) 

CALL  MULTI ( A1 , IM, IM, SAI , IM, AIN) 

CALL  MULTI ( A, IM, IM, AIN, IM, AJN) 

ERRSM=0 . ODO 
DO  5  1=1, IM 
DO  5  J=1 , IM 

ERR( I , J)=AJN( I , J) -UNI ( I , J) 

ERRSM=ERRSM+DAB  S ( ERR ( I , J) ) 

5  CONTINUE 

I F ( ERRSM . LT . 1 . OD- 7 )  GO  TO  20 
I COUNT= I COUNT + 1 
10  DO  15  1=1, IM 
DO  15  J=1 , IM 
15  A1 ( I , J ) =AIN ( I , J ) 

GO  TO  2 
20  RETURN 
END 

SUBROUTINE  INVERS (B, N, A) 

DOUBLE  PRECISION  A(N, N) , B(N, N) , C( 18 , 18 ) 
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DOUBLE  PRECISION  AMAX, TEMP, PIVOT 
DIMENSION  INDEX (4, 2 ) 

I F ( N . GT .  40 )  WRITER  101) 

101  FORMAT (2 OX, 'MATRIX  INVERSION  IS  LIMITED  TO  A 
+40  X  40  MATRIX' ) 

IF(N. GT . 40 )  GO  TO  134 

103  FORMAT (/, '  MATRIX  INVERSION'//) 

104  FORMAT('  THE  MATRIX  (XT*X)',/) 

106  FORMAT ('  ROW',  13 , IX, 1P4E16 . 7/( 8X, 1P4E16 . 7 ) ) 

128  FORMAT ( '  THE  INVERSE  OF  MATRIX  (XT*X)',/) 

131  FORMAT('  THE  UNIT  MATRIX') 

133  FORMAT('  ZERO  PIVOT') 

DO  90  1=1, N 
DO  90  J=1 ,  N 
A( I , J)=B( I , J) 

90  CONTINUE 

WRITE ( 6, 103 ) 

WRITE ( 6, 104) 

DO  105  1=1, N 

105  WRITE (6, 106)  I , (A( I, J) , J=1,N) 

DO  107  1=1 ,N 

DO  107  J=1 ,  N 

107  B  ( I ,  J )  =A  ( I ,  J ) 

DO  108  1=1,  N 

108  INDEX ( I , 1 )=0 
11=0 

109  AMAX=- 1 . 0D0 
DO  110  1=1, N 

IF( INDEX ( 1,1))  110,111,110 

111  DO  112  J=1,N 

IF( INDEX( J, 1 ) )  112,113,112 

113  TEMP=DABS ( A ( I , J ) ) 

I F ( TEMP- AMAX )  112,112,114 

114  IROW=I 
I COL= J 
AMAX=TEMP 

112  CONTINUE 

110  CONTINUE 

IF(AMAX)  225,115,116 
116  INDEX ( I COL , 1 ) = I ROW 

IF( IROW-ICOL)  119,118,119 

119  DO  120  J=1 ,  N 
TEMP=A ( I ROW , J ) 

A( IROW, J)=A( ICOL, J) 

120  A( ICOL, J)=TEMP 
11=11+1 

INDEX( 11,2 )=IC0L 
118  PIVOT=A( ICOL, ICOL) 

A( ICOL, ICOL)=l . 0D0 
PIVOT=l . ODO/PIVOT 
DO  121  J=1 , N 

121  A( ICOL, J)=A( ICOL, J ) *PIVOT 


oooooooooo 


DO  122  1=1,  N 
I F ( I - I COL )  123,122,123 

123  TEMP=A( I , I COL) 

A ( I , I COL ) =0 . 0D0 
DO  124  J=1,N 

124  A( I , J )=A( I , J ) -A( I COL, J ) *TEMP 
122  CONTINUE 

GO  TO  109 

125  I COL= I NDEX (11,2) 

IROW=INDEX( ICOL, 1 ) 

DO  126  1=1, N 
TEMP=A( I, I ROW) 

A(I, IROW)=A(I, ICOL) 

126  A ( I , I COL ) =TEMP 
11=11-1 

225  IF(II)  125,127,125 

127  WRITE( 6, 128 ) 

DO  129  1=1, N 

129  WRITE (6, 106)  I , ( A( I , J) , J=1 , N) 

DO  130  1=1, N 

DO  130  J=1,N 
C(I, J)=0.0D0 
DO  130  K=1,N 

130  C(I,J)=C(I,J)+B(I,K)*A(K,J) 

WRITE (6,131) 

DO  132  1=1, N 

132  WRITE>(6, 106)  I ,  (C(  I ,  J) ,  J=1,N) 

GO  TO  134 
115  WRITE ( 6 , 133 ) 

134  RETURN 
END 

*************************************** 

*  PROGRAM  "  ENER2"  * 

*************************************** 

THIS  PROGRAM  COMPUTES  ENERGY  DISSIPATED  IN  A 
2ND  ORDER  SYSTEM  ,  A0,A1  COEFF  GENERATED 
IN  IDENTIFICATION  PROGRAM 
DT  DELTA  T 
N  NO  OF  PTS 

DIMENSION  F( 1024) ,XINT( 1024) 

REAL  K1,K2,L1,L2 

CALL  OPSYS ( ' ALLOC ' , ' RH1 ' , 55 ) 

CALL  OPSYS ( ' ALLOC ' , ' ID1 ' , 57 ) 

READ ( 5 , * ) N 

READ(5, * )DT, AO, Al, AMS 
READ ( 55 ,90)(F(I),I=1,N) 

XINT( 1 )=0 . 0 
NM=N- 1 
Cl=- Al/AMS 
C2=-A0/AMS 


ooonooooooo 
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C3= 1.0/AMS 
Y0=0 . 0 
20=0 . 0 

DO  999  1=1 , NM 
T=DT* ( I — 1 ) 

K1=DT*Z0 

L1=DT* (C1*Z0+C2*Y0+C3*F( I ) ) 

K2=DT* ( Z0+L1 ) 

L2=DT* (Cl* (Z0+L1 ) +C2* ( Y0+K1 ) +C3*F( 1+1 ) ) 
Yl=Y0+0 . 5* (K1+K2 ) 

Zl=Z0+0. 5* (L1+L2 ) 

IP=I+1 

XINT( IP)=Z1**2 

Y0=Y1 

Z0=Z1 

WRITE ( 57,  92 )T, YO 
999  CONTINUE 

CALL  S IMP ( XINT , ENED , DT , N ) 

ENED= A 1 *  ENED + 0 . 5*A0*Y1**2 
WRITE (6, 12) ENED 
12  FORMAT (5X,7HENED  =  ,E12. 4) 

90  FORMAT (E12. 4) 

91  FORMAT (1E12. 4) 

92  FORMAT (2E12. 4) 

STOP 

END 

SUBROUTINE  SIMP(X,E,D,N) 

DIMENSION  X( 1024) 

N2=N/2-l 
E=X ( 1 ) 

DO  1  1=1, N2 
IA=2*I 
IB=IA+1 
E=E+4 . 0*X( IA) 

E=E+2 . 0*X( IB) 

1  CONTINUE 
N3=N- 1 

E=E+0 . 5*X(N3 ) +1 . 5*X(N) 

E=D*E/3 . 0 

RETURN 

END 


*************************************** 

*  PROGRAM  "  ENER3 "  * 

*************************************** 

THIS  PROGRAM  USES  THE  INCREMENTAL  EQUATIONS 
TO  COMPUTE  THE  ENERGY  DISSIPATED 
IN  A  3RD  ORDER  SYSTEM 
N  NO  OF  PTS 
DT  DELTA  T 

AO , A1 , A2  COEFF  FROM  ID  PROGRAM 


o  u  o  o 
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AO  COMPARE  TO  STIFFNESS 
ENED  TOTAL  ENERGY  DISSIPATED 
ENED1  SPRING  ENERGY  DISSIPATED 

CALL  OPSYS ( ' ALLOC ' , ' ID2 ',21) 

DIMENSION  F( 1024) ,FD( 1024) ,XINT(1024) 
DIMENSION  ESP ( 1024) , XINT1 ( 1024 ) 

READ ( 5 , * )N 

READ(5, * )DT,A0, Al, A2, AMS 
CALL  OPSYS ( ' ALLOC  1 ,  ' RH1 ' , 55 ) 

111  READ( 55,91)(F(I), FD( I ) , 1=1 , N) 

C  CONVERT 

A0=A0/(A2*AMS) 

A1=A1/ ( A2  * AMS ) 

B1=1.0/(A2*AMS) 

B2=l. 0/AMS 
A2=l . 0/A2 
C  CONSTANTS 

R0=2  4 . 0/DT*  *  3 
Rl=4.0 
R2=12 . 0/DT 
R3=24.0/DT**2 
SO=12 . 0/DT* *2 
S1=DT 
S2=6.0 
S3=12 . 0/DT 
U0=4 . O/DT 
U1=DT** 2/6.0 
U2=DT 
113=4 . 0 

Ql=24. 0/DT**3+12 . 0*A2/DT**2+4 . 0*A1/DT+A0 

Q2=4.0+A2*DT  A1*DT* *2/6.0 

Q3=12 . 0/DT +6 . 0*A2+Al*DT 

Q4=24 . 0/DT*  *2  + 12 . 0*A2/DT+4 . 0* Al 

Y0=0 . 0 

Y1=0 . 0 

Y2=0 . 0 

Y3=B1*F( 1 ) +B2*FD( 1 ) 

XINT( 1 )=0 . 0 
ESP ( 1 )=0 . 0 
NM=N-1 

DO  999  1=1, NM 
T=DT* ( I — 1 ) 

DF=B 1*(F(I+1)-F(I) )+B2*(FD(I+l)-FD(I) ) 

DZ0=(DF+Q2*Y3+Q3*Y2+Q4*Y1)/Q1 

DZ1=U0*DZ0-U1*Y3 -U2  * Y2 -U3  *Y1 

DZ2=S0*DZ0-S1*Y3-S2*Y2-S3*Y1 

Y0=Y0+DZ0 

Y1=Y1+DZ1 

Y2=Y2+DZ2 

Y3=B1*F(I+1) +B2*FD( 1+1 )-A2*Y2-Al*Yl-A0*Y0 
IP=I+1 


noon  noonooo 
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XINT ( IP )=Y1* ( F ( IP ) -AMS* Y2 ) 

WRITE (21,91)1, YO 
999  CONTINUE 

CALL  SIMP(XINT,ENED,DT,N) 

WRITE ( 6 , 12 )ENED 
12  FORMAT ( 5X ,  7HENED  =  ,E12.4) 

90  FORMAT (E12. 5) 

91  FORMAT (2E12. 4) 

STOP 

END 

SUBROUTINE  SIMP { X( E , D, N ) 

DIMENSION  X( 1024 ) 

N2=N/2-l 

E=X(1) 

DO  1  1=1 ,N2 
IA=2*I 
IB=IA+1 
E=E+4 . 0*X( IA ) 

E=E+2 . 0*X( IB ) 

1  CONTINUE 
N3=N-1 

E=E+0.5*X(N3)+1.5*X(N) 

E=D*E/3 . 0 

RETURN 

END 

****************************************** 

*  PROGRAM  "  FREQ ID"  * 

****************************★**-*********** 


FREQUENCY  DOMAIN  SYSTEM  IDENTIFICATION 

REAL  DD ( 1 02 4 ) , F ( 1 02  4 ) , W ( 1 02  4 ) , COE ( 4 ) , QT ( 2  0 ) 
DIMENSION  S ( 1024 ) , INDEX( 1024 ), QQ( 1024 ), X2 ( 22 , 1 ) 
DIMENSION  XI ( 22 , 2  )  , XT1 ( 2 , 22 ) , XP ( 2 , 2 ) , XPI ( 2 , 2 ) 
DIMENSION  XTT ( 2 , 2  2  )  ,  B  (  2  , 1) 

COMPLEX  FT( 1024) ,DDF( 1024) 

CALL  OPSYS ( ' ALLOC ' , ' FD1 ' , 10 ) 

CALL  OPSYS ( ' ALLOC ' , ' ID1 ' , 8 ) 

COMMON  /ARRAY/W , QQ , S , INDEX 
EXTERNAL  SPLINE 
EXTERNAL  QF 
EXTERNAL  MATINV 

WHEN  INDE=1  SHOWS  THE  PROGRAM  WILL  GO  TO 
THE  THIRD  ORDER  SYSTEM 

INDE=1 

ID=0 

KDEX=0 

MP=2 

M=10 


oooooo  non  non  ooo  ooo 


N=1024 
NT=N/2 
NTP=NT+1 
DT=0 . 05 
C2=l . 0 

PI=3. 1415926535 

NS=N/2 

CA=P I/NS 

GET  THE  INPUT  FROM  THE  FILE  AND  FFT  THE  INPUT 

READ(10,15)  ( F( I ) , DD ( I ) , 1=1 , N) 

15  FORMAT (2E 12. 4) 

I F ( KDEX . NE . 1 )  GO  TO  13 

SMOOTH  THE  FUNCTION 

DO  9  1=1, NS 
IM=N- 1+1 

CMULT=0 . 5* ( 1 . O-COS (  (  I-1)*CA) ) 

F(IM)=CMULT*F(IM) 

DD ( IM ) =CMULT*DD ( IM ) 

9  CONTINUE 
13  TT=N*DT 

DO  THE  FOURIER  TRANSFORM 
DO  20  1=1, N 

FT ( I ) =CMPLX (F(I),0.0) *TT 
DDF ( I )=CMPLX(DD( I ) , 0 . 0 ) *TT 
20  CONTINUE 

CALL  FFT1 ( FT , M , N , - 1 . ) 

CALL  FFT1 (DDF, M, N, -1.0) 

CALCULATE  THE  DF  AND  FREQUENCY 

DF=2 . 0*PI/TT 
DO  24  1=1, N 
W(I)=(I-1) *DF 
24  CONTINUE 

IF(ID.EQ.l)  GO  TO  33 
GO  TO  45 

IDENTIFICATION  PROCESS 
SECOND  ORDER  SYSTEM 

CHANGE  F.T.  OF  ACCELERATION  TO  F.T.  OF  DISPLACEMENT 

33  DO  55  1=2, N 

DDF ( I )=-DDF (I)/(W(I)**2) 

55  CONTINUE 
45  DO  50  1=40,60 


ono  ooooooooo  ooo  ooo  ooo 
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FF=CABS ( FT( I ) ) 

ZF=CABS { DDF ( I ) ) 

QQ(I)=FF**2/(ZF**2) 

WRITE ( 8, 88 ) W( I ) , QQ( I ) 

50  CONTINUE 
GO  TO  200 

FIND  THE  MINIMUN  OF  THE  QQ(I)  AND  IT'S  FREQUENCY 

SMAL=QQ(30) 

KN=30 

DO  60  1=31,70 
I F ( SMAL . LE . QQ ( I ) )  GO  TO  60 
SMAL=QQ ( I ) 

QMIN=QQ( I ) 

KN=I 

60  CONTINUE 
WM=KN*DF 

CALCULATE  THE  STIFFNESS 
CO=C2*WM**2 

WRITE ( 6, * )  KN,  WM, CO, QMIN 

USE  FIT  TO  FIND  THE  REAL  MINIMUM 

AS=0.9*WM 
AF=1.1*WM 
I S=AS/DF 
I F=AF/DF 
IJ=IS-1 
ITT=(IF-IS)+1 
DO  67  1=1, ITT 
I K= I  + 1 J 
QT(I)=QQ(IK) 

67  CONTINUE 

POLYNOMIAL  FITTING  BY  LEAST  SQUARE  METHOD 

QT  : INPUT  VALUES 

ITT: TOTAL  PTS  OF  QT 

AS -.START  FREQUENCY 

AF : FIANL  FREQUENCY 

COE( I ) : COEFFICIENTS  OF  POLYNOMIAL 

CALL  FIT1 (QT, 3 , ITT, AS , AF, COE ) 

WM=-COE( 2 )/( 2 . 0*COE ( 3 ) ) 

CO=C2*WM**2 

WRITE ( 6 , * ) ITT, WM, CO 

SET  THE  RANGE  OF  FREQUENCY  AND  FIT  THE 
DATA  BY  CUBIC  SPLINE 


QA=0 . 9 
QB=1 • 1 
WA=QA*WM 
WB=QB*WM 

IF( INDE.EQ. 1)  GO  TO  80 
CALL  SPCOEF(NT, W, QQ,  S, INDEX) 

DO  70  1=1, NT 
QQ(I)=-QQ(I) 

70  CONTINUE 

DO  THE  INTEGRATION  AND  CALCULATE  THE  DAMPING 

CALL  SIMP (QF , WA, WB, 1 . OE-5 , ANS, ERROR, AREA, I FLAG) 
QIN=5.*ANS/(WM**7) 

Q3=(QB**3)-(QA**3) 

Q5=(QB**5)-( Q A  *  *  5 ) 

Q7=(QB**7)-(QA**7) 

QP=(-5./3 . )*Q3+2. *Q5-(5./7. )*Q7 
CI=( (WM**2)/Q5)*(QIN+(C2**2)*QP) 

CC=SQRT(CI) 

WRITE  (6,*)  CO, CC, ANS, I FLAG 
GO  TO  200 

IDENTIFICATION  FOR  THIRD  ORDER  SYSTEM 

CALCULATE  THE  FREQUENCY  BAND  FOR  IDENTIFICATION 

80  IA=WA/DF 
IB=WB/DF 
I  I=IA-1 
NF=( IB-IA) +1 
WRITE ( 6, * )NF, IA, IB 

FORMING  THE  IDETI FI CATION  MATRIX 

DO  100  I=IA, IB 

J=I-II 

EU=W( I ) **2 

EV=-2 . *C2* (W( I ) **4) 

EE=WM**2- (W(I)**2) 

EW=(C2**2 ) * (EE* *2 ) 

QQ ( I ) =COE ( 1 ) +COE ( 2 ) *W ( I ) +COE (3)*(W(I)**2) 

XI (J, 1 )=EU 
XI (J, 2 )=EV 
X2 ( J , 1 ) =QQ ( I ) -EW 
100  CONTINUE 

DO  105  1=1,2 
DO  105  J=1 , 2 
XP(I, J)=0.0 
DO  105  K=1 , NF 

XP(I, J)=XP(I, J)+X1(K, I ) *X1 (K, J) 


105  CONTINUE 

DO  104  1=1,2 
DO  104  K=1,NF 
XT1 ( I , K)=X1 (K, I ) 

104  CONTINUE 

CALL  MATINV(XP,MP,XPI) 

DO  106  1=1,2 

WRITE ( 6,15)  (XPI(I, J) , J=l,2) 

106  CONTINUE 

DO  110  1=1,2 
DO  110  J=1 , NF 
XTT ( I , J ) =0 . 0 
DO  110  K=l, 2 

110  XTT( I , J)=XTT( I , J) +XPI ( I , K) *XT1(K, J) 

DO  118  1=1,2 

118  WRITE( 6, 37 ) (XTT( I , J) , J=l, NF) 

DO  115  1=1,2 

B ( I , 1 ) =0 . 0 

DO  115  K=1,NF 

115  B ( I , 1 ) =B (1,1) +XTT ( I , K ) *X2 ( K , 1 ) 

DO  119  1=1,2 

119  WRITE(6, * )  B( I , 1 ) 

CK=SQRT(B(1, 1) ) 

CL=B(2,1)/CK 
WRITE ( 6, * ) CK, CL 

35  FORMAT ( 2 4H  REAL  PART  I MAG INAL  PART) 

37  FORMAT (6E 12. 4) 

88  FORMAT (2E12. 4) 

200  STOP 
END 

FUNCTION  QF  IS  A  CONTINUOUS  FUNCTION. 

THE  DATA  QQ(I)  WAS  ARRANGED  SO  THAT  THE  FUNCTION 
COULD  BE  CALLED  IN  ANY  TIME. 

FUNCTION  QF(WW) 

DIMENSION  W(1024) ,QQ( 1024), S( 1024), INDEX( 1024) 

COMMON  /ARRAY/W,QQ,S, INDEX 

N=512 

X=WW 

QT=SPLINE(N,W,QQ,S, INDEX, X) 

QF=QT* (X**2 ) 

RETURN 

END 
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*  PROGRAM  NAME  "PUR"  * 

****************************************** 

THIS  PROGRAM  IDENTIFIES  THE  PARAMETERS  OF  A 
TIME  VARYING  SECOND-ORDER  LINEAR  MODEL  BY 
USING  THE  PERTURBATION  METHOD  AND  ITERATIVE 
NEWTON-  RAPHSON  PROCEDURE.  IT  MAY  ALSO  USE 
THE  POLYNOMIAL  FITTING  APPROACH. 

COMPLEX  FF( 1024) , FD( 1024) 

COMPLEX  ZZ( 1024) ,HD,Z(1024) /HABS( 512) 
DIMENSION  DD( 1024) , F( 1024 ) ,HMOD( 512 ) , HZ( 512 ) 
DIMENSION  EK(3),C0E(4),S(9), INDEX( 9 ) , AS( 9 ) 
COMMON/ARRAY/FF , FD 
CALL  OPSYS ( ' ALLOC ' , ' FD1 ' , 10 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH1 '  ,  7 ) 

CALL  OPSYS ( ' ALLOC ' , ' NFL ' , 4 ) 

EXTERNAL  MATINV 

SM  :  SYSTEM  MASS 

ALPHA  :  DAMPING  COEFFICIENT 

BETA  :  STIFFNESS  COEFFICIENT 

ZETA  :  DAMPING  RATIO 

Cl  :  INITIAL  GUESS  'DAMPING' 

C2  :  INITIAL  GUESS  'ALPHA*C1' 

C3  :  INITIAL  GUESS  'STIFFNESS' 

C4  :  INITIAL  GUESS  ' BETA*C3 ’ 

ESLON  :  ACCURACY  MEASURE  (FOR  ESLON1-ESLON4*) 
W1-W4  :  ACCURACY  MEASURE  FOR  NEWTON  METHOD 
DA1-DA4  :  INCREMENT  VALUES 

SM=1 . 0 
Cl=l .  5 
C2=0 . 02 
C3=37 . 5 
C4=-0 . 01 
ZETA=0 . 1 
ESLON1=0 . 02 
ESLON2=0 . 02 
ESLON3=0 . 02 
ESLON4=0 . 02 
WN=SQRT ( C3/SM ) 

DA1=0 . 01 
DA2 =0.001 
DA3=1 . 0 
DA4"0.01 
Wl=0.01 
W2=0 . 01 
W3=0 . 01 
W4=0 . 01 

PI=3. 1415926535 
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DT=0 . 05 
M=10 
N=2**M 
NT=N/8 
TT=DT*N 
DW=2 . *PI/TT 
NS=N/2 
CA=PI/NS 
C 

C  READ  IN  THE  TIME  DOMAIN  RESPONSE 
C 

READ( 10,80)  (F(I)/DD(I),I=1/N) 

C 

C  WINDOW  THE  DATA 

C 

DO  3  1=1, NS 
IM=N- 1+1 

CMULT=0 . 5* ( 1 . O-COS ( ( 1-1 ) *CA) ) 
F(IM)=CMULT*F(IM) 

DD ( IM ) =CMULT*DD ( IM ) 

3  CONTINUE 
C 

C  FFT  THE  INPUT  AND  RESPONSE 

C 

DO  8  1=1, N 

FF ( I ) =CMPLX ( F ( I ) , 0 . 0 ) *TT 
ZZ ( I ) =CMPLX ( DD ( I ) , 0 . 0 ) *TT 

8  CONTINUE 

CALL  FFT1 (FF,M,N, -1.0) 

CALL  FFT1(ZZ,M,N, -1.0) 

C 

C  FIND  THE  MOD ( H ) 

C 

DO  9  1=1, NT 
FA=CABS ( FF ( I ) ) 

ZA=CABS ( ZZ ( I ) ) 

HMOD( I )=ZA/FA 

9  CONTINUE 
C 

C  CALCULATE  THE  DERIVATIVE  OF  FF 

C 

DO  12  1=1, N 
TC= ( I - 1 ) *DT 
F( I )=F ( I ) *TC 

FD( I )=CMPLX(F( I ) , 0 . 0 ) *TT 

12  CONTINUE 

CALL  FFT1(FD,M,N, -1 . ) 

DO  13  1=1, N 

13  FD ( I ) =FD ( I ) *CMPLX ( 0 . 0 , - 1 . ) 

C 

C  FIND  THE  NUMBER  OF  PTS  NEEDED  IN  THE  BAND 
C 
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NC= (1.0*  ZETA* WN ) /DW 
NW=WN/DW 
NS=NW-NC 
NF=NW+NC 
AA1=C1 
AA2=C2 
AA3=C3 
AA4=C4 
26  L=1 

28  WRITE (6,89)01,02,03, C4 
IC0UNT=2 

I STEA= I COUNT 

AC1=C1 

AC2=C2 

AC3=C3 

AC4=C4 

CALL  SUBROUTINE  PURT  TO  FIND  THE  MODULUS 
OF  H  FROM  THE  MEASURED  INPUT 

29  WN=SQRT ( C3/SM ) 

NC=( 1 . 0*ZETA*WN)/DW 
NW=WN/DW 
NS=NW-NC 
NF=NW+NC 

CALL  PURT ( DW , NT , Cl , C2 , C3 , C4 , HZ ) 

FIND  E( ERROR  MEASURE) 

EE=0 . 0 

DO  30  I=NS,NF 
E=HMOD ( I ) -HZ ( I ) 

EE=EE+(E**2) 

30  CONTINUE 

EK ( I COUNT ) =EE 
I F ( I COUNT . EQ . 3 )  GO  TO  50 
GO  TO  (31, 32, 33, 34), L 

31  IF( ICOUNT . EQ . 1 )  GO  TO  40 
C1=C1-DA1 

I COUNT= I COUNT- 1 
GO  TO  29 

40  C1=AC1+DA1 

I COUNT= I STEA+ 1 
GO  TO  29 

32  IF( ICOUNT. EQ. 1)  GO  TO  41 
C2=C2-DA2 

I COUNT= I COUNT - 1 
GO  TO  29 

41  C2=AC2+DA2 
ICOUNT=ISTEA+l 
GO  TO  29 

33  IF( ICOUNT. EQ. 1)  GO  TO  42 
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C3=C3-DA3 
I COUNT= I COUNT- 1 
GO  TO  29 

42  C3=AC3+DA3 
IC0UNT=ISTEA+1 
GO  TO  29 

34  I F ( I COUNT . EQ . 1 )  GO  TO  43 
C4=C4-DA4 
I COUNT= I COUNT - 1 
GO  TO  29 

43  C4=AC4+DA4 
ICOUNT=ISTEA+l 
GO  TO  29 

CALCULATE  THE  FIRST  AND  SECOND  DERIVETIVES 

50  GO  T0(60, 61, 62 , 63 ) , L 

60  TD=DA1 
AV=AC1 
W=W1 

GO  TO  70 

61  TD=DA2 
AV=AC2 
W=W2 

GO  TO  70 

62  TD=DA3 
AV=AC3 
W=W3 

GO  TO  70 

63  TD=DA4 
AV=AC4 
W=W4 

GO  TO  70 

70  DEK=(EK(3 ) -EK(  1 ) ) / ( 2  - 0*TD) 

DDEK=(EK(3)-2.*EK(2)+EK(1) )/(TD**2) 

DO  THE  NEWTON  RAPHSON  METHOD 

ACC=AV- ( DEK/DDEK ) 

WE=ABS ( ( ACC- AV ) /ACC ) 

WD=WE-W 

IF(WD.LE.O.O)  GO  TO  87 
99  GO  T0( 101 , 102 , 103 , 104) , L 

101  C1=ACC 
GO  TO  28 

102  C2-ACC 
GO  TO  28 

103  C3=ACC 
GO  TO  28 

104  C4=ACC 
GO  TO  28 

87  GO  TO(92,93,94,95)  ,L 
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92  C1=ACC 
GO  TO  85 

93  C2=ACC 
GO  TO  85 

94  C3=ACC 
GO  TO  85 

95  C4=ACC 
85  L=L+1 

IF(L.EQ.5)  GO  TO  100 
GO  TO  28 

COMPARE  THE  PARAMETERS 

100  CD1=ABS( (AA1-C1)/C1) 

CD2=ABS ( ( AA2-C2 ) /C2 ) 

CD3=ABS ( ( AA3-C3 ) /C3 ) 

CD4=ABS ( ( AA4-C4 ) /C4 ) 

IF(CD1 . GT. ESLON1 )  GO  TO  81 
I F ( CD2 . GT . ESLON2 )  GO  TO  81 
I F ( CD3 . GT . ESLON3 )  GO  TO  81 
I F ( CD3 . GT . ESL0N4 )  GO  TO  81 
WRITE(6, * )  01,02,03,04 
GO  TO  88 

81  AA1=C1 
AA2=C2 
AA3=C3 
AA4=C4 
GO  TO  26 

80  FORMAT (2E12. 4) 

82  FORMAT (3E13. 4) 

89  FORMAT (4E 12. 4) 

150  FORMAT (2E16. 6) 

88  STOP 
END 

PERTURBATION  METHOD 

SUBROUTINE  PURT(DW, NT, Cl , 02 , 03 , 04, HA) 
COMPLEX  FF( 1024) , FD( 1024) ,HABS(512 ) , Z(512 ) 
COMPLEX  HD, H, HDEV, HI , H2 , H3 , H4 
DIMENSION  HA (NT) 

COMMON/ARRAY/FF, FD 
SM=1 . 0 

DO  14  1=1, NT 
W= ( I - 1 ) *DW 

HD=CMPLX( (C3-SM* (W**2 ) ) , W*C1 ) 

H=1 . O/HD 

HD E V=CMP LX ( ( 2 . * SM * W ) , - C 1 ) / ( HD *  * 2 ) 

H1=H*FF ( I ) 

H2=C2*H* (Hl+W* (HDEV*FF( I ) +H*FD( I ) ) ) 
H3=CMPLX(0 . 0 , -04 ) *H* (HDEV*FF ( I ) +H*FD( I ) ) 

Z( I)=H1+H2+H3 
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H4=HDEV+H* ( FD ( I ) /FF ( I ) ) 

HABS ( I ) =H* ( 1 . 0+C2 * ( H+W*H4 ) -CMP LX ( 0 . 0 , C4 ) *H4 ) 
HA( I )=CABS(HABS( I ) ) 

14  CONTINUE 
GO  TO  10 
10  RETURN 
END 

****************************************** 

*  PROGRAM  "MDFID"  * 

****************************************** 


IDENTIFICATION  OF  MDF  SYSTEM 

DIMENSION  F2 ( 1024) , Yl( 1024) , Y2 ( 1024) , W( 50) , CT( 6 ) 
COMPLEX  Y1FT( 1024) , Y2FT( 1024) , F2FT( 1024) , AINV(2 , 2 ) 
COMPLEX  B(2,2),H(2,2), H12 , H22 ,  A(  2 , 2  ) ,  B1 ,  B2 ,  CD 
COMPLEX  AA ( 2 , 2 ) , BB ( 2 , 2  )  , AAIN( 2 ,2), Bll , B22 
DIMENSION  EPS1 ( 50 ) , EPS2 ( 50 ) , Y1T ( 100 ) , Y2T ( 100 ) , CC ( 6 ) 
DIMENSION  C( 6) , DC( 6) , Yll ( 100 ) , Y22 ( 100) , ERROR ( 100) 
DIMENSION  CCOUNT ( 6 ) , E2 ( 20 ) 

REAL  Ml , M2 , MC , COE ( 5 ) 

CALL  OPSYS ( ' ALLOC ' , ' RH3 ' , 14) 

CALL  OPSYS ( ' ALLOC ' , ' MDF1 ' , 10 ) 

C  CALL  OPSYS ( ' ALLOC ' , ' MODI ' , 13 ) 

EXTERNAL  INVERS  . 

C 

C  SET  PARAMETERS 

C 

C  DT:  TIME  STEP 

C  DW: DELTA  W 

C  NT: NUMBER  OF  TIME  STEP 

C  Ml: MASS  FOR  1ST  DEGREE  OF  FREEDOM 

C  M2: MASS  FOR  2ND  DEGREE  OF  FREEDOM 

C  NWA: NUMBER  OF  FFT  FREQ.  WHERE  ANALYSIS  STARTS 

C  NID: NUMBER  OF  PARAMETERS  IDENTIFIED 

C  NW: NUMBER  OF  FREQUENCIES  WHERE  ANALYSIS  DONE 

C  NCUV: NUMBER  OF  PTS  TO  DEFINE  THE  CURVE 

C  W1 : 1ST  FREQUENCY  WHERE  ANALYSIS  DONE 

C  NCY: NUMBER  OF  ITERATION  WILL  BE  PERFORMED 

C  C( J) , J=l, 6  : SYSTEM  PARAMETERS 

C  DC( J) , J— 1, 6  : INCREMENT  FOR  C(J) 

C  I ORDER:  3=THIRD-0RDER;  2=SECOND-ORDER 

C 

NCOL=2 

IORDER=3 

ISMOTH-1 

DT=0 . 05 

NT=1024 

N3=1024/2 

NCY=20 

NID=6 
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NCUV=7 
M=10 
Ml=l . 0 
M2=l . 0 

INITIAL  GUESS 

C( 1 )=30 . 0 
C(2 )=1 . 5 
C(3 )=0 . 01 
C(4)=25 . 0 
C  ( 5 ) =2 . 0 
C( 6 )=0 . 01 
DC( 1 )=1 . 0 
DC(2)=0.05 
DC( 3 )=0 . 001 
DC(4)=1 . 0 
DC( 5 )=0 . 05 
DC( 6 )=0 . 001 
ACUR=0 . 02 
PI=3. 1425926 
CA— P I /NS 
TT=DT*NT 
DW=2.0*PI/TT 
NW=30 
NWA=30 
Wl=NWA*DW 
DO  111  1=1, NW 

111  W( I )=W1+ ( 1-1 ) *DW 
W2=W(NW) 

DO  112  1=1, NID 

112  CT ( I ) =C ( I ) 

MC=M1+M2 

READ  INPUT  AND  REPONSE  THEN  FFT 

READ( 10,12)  ( Y1 ( I ) , Y2 ( I ) , 1=1 , NT) 

READ( 14,16)  ( F2 ( I ) , 1=1 , NT) 

COMPUTE  THE  RELATIVE  DISPLACEMENT  AND 
SMOOTH  THE  INPUT  FUNCT I ONWHEN  IF  NECESSARY 


DO  37  1=1, NT 
37  Y2 ( I ) =Y2 ( I ) -Y1 ( I ) 

IF( ISMOTH.NE. 1 )  GO  TO  39 

DO  36  1=1, NS 

IM=NT-I+1 

CMULT=0 . 5* ( 1 . O-COS ( ( 1-1 ) *CA) ) 
Y1 ( IM) =CMULT*Y1 ( IM ) 
Y2(IM)=CMULT*Y2(IM) 

F2( IM)=CMULT*F2( IM) 

36  CONTINUE 
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IF( IORDER . EQ . 2 )G0  TO  200 
B1=CMPLX( 1 . 0, W( IW) *CT( 3 ) ) 

B2=CMPLX(1.0,W(IW)*CT(6) ) 

Bll=( W( IW) ) **2*MC*B1 
B22=(W( IW) ) **2*M2*B2 

A(  1 , 1 )=CMPLX(CT( 1 ) -REAL ( Bll ) , CT ( 2 ) *W( IW) -AIMAG( Bll ) ) 
A(1,2)=-(W( IW) ) **2*M2*B1 
A(2, 1 )=-B22 

A(2 , 2 )=CMPLX(CT( 4) -REAL(B22 ) , CT( 5 ) *W( IW) -AIMAG(B22 ) ) 
CALL  MATINV(A,2, AINV) 

B( 1 , 1 )=B1 
B(  1 , 2 )=B1 

B(2, 1)=CMPLX(0. 0,0.0) 

B( 2 , 2 )=B2 

CALL  MULTI(AINV,2,2,B,NCOL,H) 

GO  TO  205 

CALCULATION  FOR  SECOND  ORDER 

200  AA( 1,1) =CMPLX( CT( 1 ) -MC* ( W ( IW) ) **2 , CT( 2 ) *W( IW) ) 

AA(2 , 2 )=CMPLX(CT( 3 ) -M2  *(W(IW))**2, CT( 4 ) *W( IW) ) 
AA(1,2)=CMPLX(-M2*(W( IW) ) **2,0.0) 

AA (2,1) = AA (1,2) 

BB (1,1) =CMPLX ( 1 . 0 , 0 . 0 ) 

BB( 1 , 2 )=CMPLX( 1 . 0 , 0 . 0 ) 

BB ( 2 , 1 ) =CMPLX( 0 .0,0.0) 

BB (2,2) =CMPLX ( 1 . 0 , 0 . 0 ) 

CALL  MATINV(AA, 2 , AAIN) 

CALL  MULTI ( AAIN , 2,2, BB , NCOL , H ) 

205  H12=H( 1,2) 

H22=H( 2 , 2 ) 

Y1T ( I W ) =CABS ( H12 * F2FT ( NWA+ I W- 1 ) ) 

Y2T( IW)=CABS(H22*F2FT(NWA+IW-1) ) 

Yll ( IW)=Y1 (NWA+ IW-1 ) 

Y22 ( I W ) =Y2 ( NWA+ I W- 1 ) 

CALCULATE  THE  ERROR  TERM 

EPS1( IW)=Y11( IW)-Y1T( IW) 

EPS2 ( IW) =Y22 ( IW) -Y2T( IW) 

96  CONTINUE 
E2 ( IX)=0 . 0 

SUMMATION  OF  ERROR 

DO  121  1=1 , NW 

121  E2 ( IX)=E2 (IX) +EPS1 ( I)**2+EPS2(I)**2 
WRITE (6,12)  CT( IDX) , E2 ( IX) 

97  CONTINUE 
WRITE (6, 25) 

FIT  THE  E2 ( IX)  IN  A  POLYNOMIAL  EQUATION  BY 
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USING  THE  LEAST  SQUARE  METHODAND  THEN  FIND 
THE  CORRESPONDING  FREQUENCY  WHERE  THE  CURVE 
HAS  A  MINIMUM 

CALL  FIT1 ( E2 , 3 , NCUV, CTI , DC ( IDX) ,COE) 

CHECK  THE  CUVATURE 

103  CT( IDX) =-C0E ( 2 )/( 2 . 0*COE ( 3 ) ) 

C( IDX)=CT( IDX) 

CCOUNT { I COUNT )=CT( IDX) 

I C0UNT= I COUNT + 1 

98  CONTINUE 

ERROR ( INDEX ) =E2 ( NCUV ) 

ACU=ABS( ERROR ( INDEX- 1 ) -ERROR ( INDEX) ) /ERROR ( INDEX) 
IF(ACU.LT. 1.0E-3)  GO  TO  100 
DO  191  1=1 , NID 

RATIO=ABS( CCOUNT ( I ) -CC( I ) )/CC( I ) 

IF (RATIO . LE . ACUR )  GO  13  191 
GO  TO  95 
191  CONTINUE 
GO  TO  100 

95  WRITE ( 6,19)(C(I),I=1,NID) 

99  CONTINUE 

100  WRITE ( 6, * ) ( C ( I ) , 1=1 , NID ) 

GO  TO  87 

12  FORMAT (2E12. 4) 

16  FORMAT (E12. 4) 

17  FORMAT (3E12. 4) 

19  FORMAT (6E12. 4) 

23  FORMAT ( 4E1 2 . 4 ) 

25  FORMAT (/) 

87  STOP 
END 

MATRIX  MULTIPLICATION 

SUBROUTINE  MULTI (TA, N, M, TB, L, TC) 

COMPLEX  TA(N,M),TB(M,L),TC(N,L) 

DO  15  1=1, N 
DO  15  J=1 , L 
TC ( I , J ) =0 . 0 
DO  15  K=1 , M 

15  TC ( I , J ) =TC ( I , J ) +TA ( I , K) *TE ( K , J ) 

RETURN 

END 
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LINEAR  RESPONSE  OF  HIGH  ORDER  MDF  SYSTEM 


MM  -HIGHEST  DERIVATIVES 

NT  ^NUMBER  OF  POINTS  ANALYZED 

DT  =TIME  STEP 

A(J,K,L)=(AJ) 

X0( J, K)=JTH  DERIVATIVES  OF  STATE  VECTOR  AT  TIME  T 
DX ( J ,  K  ) “DELTA  OF  JTH  DERIVATIVES 

XN( J, K)=JTH  DERIVATIVES  OF  STATE  VECTOR  AT  TIME  T+DT 
F(J,K)  =JTH  FORCE  AT  TIME  K 

DIMENSION  A(4, 2,2) ,X0(4,2 ) ,DX(4,2) ,XN(4,2) ,CA(2,2) 
DIMENSION  CB(2,2) ,  F(2, 1024) , DM(2,2) , AMI (2, 2) ,FT(5) 
DIMENSION  DV(2 ) , DU(2 ) , DW(2 ) , DN(2 , 2 ) , DR(2 ) , FF( 1024) 
DIMENSION  FK(5),DD(4,2),DE(4,2),DF(2),FD(1024) 
DIMENSION  ENED1 ( 1024 ) , ENED2 (1024), FDD (2,1024) 

CALL  OPSYS ( ' ALLOC' , ' RH1 ' , 10 ) 

CALL  OPSYS ( 'ALLOC' , 'MDF2' ,16) 

SPECIFY  DATA(AJ)=A(J,K,L) 

A(l, 1, 1)=29.07 

A(l, 1,2)=0.0 

A(l,2,l)=0.0 

A(l,2,2)=32.3 

A(2, 1, 1)=2. 11 

A(2,l,2)=0.0 

A(2,2,l)=0.0 

A(2,2,2)=1.89 

A(3,l,l)=l. 0+1.0 

A(3, 1,2)=1.0 

A(3,2,l)=1.0 

A(3,2,2)— 1.0 

A(4, 1, 1)=0.0082*A(3, 1,1) 

A(4, 1, 2)=0 . 0082*A( 3 ,1,2) 

A(4,2,1)=0.0048*A(3,2, 1) 

A(4, 2, 2 )=0 . 0048*A(3 , 2,2) 

M=4  MEANS  THIRD-ORDER;  M=3  MEANS  SECOND-ORDER 


M-4 
MM=M-1 
DT=0 . 05 
NT=1024 
ENED1 ( 1 )=0 . 0 
ENED2 ( 1 )=0 . 0 
AM1=1 . 0 
AM2=1 . 0 
DO  14  1=1, M 
DO  14  J=1 , 2 
X0(I,J)=0.0 


14  CONTINUE 


READ  THE  INPUT 

READ(10, 18) (FF( I ) , FD( I ) , 1=1, NT) 
18  FORMAT (2E12. 4) 

DO  20  1=1, NT 
F ( 2 , I ) =FF ( I ) 

F(1,I)=FF(I) 

20  CONTINUE 

IF(M.NE.4)  GO  TO  91 

DO  93  1=1, NT 

FDD(1, I )=FD( I ) *A(4, 1,2) 

FDD(2 , I )=FD( I ) *A(4,2 , 2 ) 

93  CONTINUE 


FIND  A(M)  AND  ITS  INVERSE 

91  DO  121  1=1,2 
DO  121  J=l, 2 
121  DM(I, J)=A(M, I, J) 

CALL  MATINV(DM,2,AMI) 

FACTOR  COMPUTATION 


FT(M)=1 . 0 
DO  131  J=2,M 
JJ=M+1-J 

131  FT( JJ)»FT( JJ+1)*FL0AT( J) 

CALCULATE  THE  CONSTANT  VALUES  OF  EQ  13 

DO  141  K=l, 2 
DO  141  L*l,2 
141  CA(K,L)=0.0 
DO  151  J=1,MM 
AC=(DT**(M-J) )/FT( J) 

DO  150  K=l,2 
DO  150  L=l,2 

150  CA(K,L)=CA(K,L)+AC*A(J,K,L) 

151  CONTINUE 

DO  161  K=l, 2 
DO  161  L=1 , 2 
161  DM(K,L)=CA(K,L) 

CALL  MULT( AMI ,DM, CA, 2 ) 

DO  171  K=l, 2 
171  CA(K,K)*CA(K,K)+1.0 
CALL  MATINV(CA,2,DM) 

CALL  MULT (DM, AMI ,CA, 2 ) 

FACTOR  COMPUTATION 
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FK( 1 )=1 . 0 

DO  173  J=2, 5 

FK(J)*FK(J-l)*FLOAT(J-l) 

173  CONTINUE 

START  THE  COMPUTATION  LOOP 

WRITE(16,26)T,X0(1,1),X0(1,2) 
DO  999  INDEXS2 , NT  . 

T=( INDEX- 1)*DT 


SOLVE  EQ  13 

(CALCULATE  THE  DIFFERENCES  OF  MTH  DERIVATIVES) 

DO  181  J=1 , 2 
181  DW(J)=0.0 

DO  191  J=1,MM 

LIM=M-J 

DO  183  L=1 , 2 

183  DV(L)*0.0 

DO  185  K*1 , LIM 
AC=(DT**K)/FK(K+1) 

JK=J+K 
DO  184  L=l,2 

184  DV(L)=DV(L)+AC*XO(JK,L) 

185  CONTINUE 

DO  186  K=l,2 
DO  186  L*l,2 

186  DN(K,L)SA( J,K,L) 

CALL  MULV(DN,DV,DU/2) 

DO  187  1*1,2 

187  DW( I )*DU( I )+DW( I ) 

191  CONTINUE 

DO  201  1*1,2 

201  DF(I)*F(I, INDEX) -F( I, INDEX-1 )+ 

+(FDD( I, INDEX)- FDD ( I , INDEX- 1 ) ) 

DO  211  1*1,2 
211  DV( I )*DF( I )-DW( I ) 

CALL  MULV(CA,DV,DW,2) 

DO  221  1*1,2 
221  DX(M, I )*DW( I ) 

SOLVE  EQ  7 

(CALCULATE  THE  DIFFERENCES  OF  JTH  DERIVATIVES) 

DO  261  J*1,MM 
DO  230  1*1,2 
230  DV( I )*0.0 
LIM*M-J 

DO  231  K*l, LIM 
AC*(DT**K)/( FK(K+1 ) ) 


ooo  oooooo  oooo  ono 


DO  228  L*l,2 

228  DV(L)*DV(L)+AC*XO(JK,L) 

231  CONTINUE 

AC=(DT**(M-J) )/FT( J) 

DO  241  1*1,2 

241  DX( J, I )*DV( I )+AC*DX(M, I ) 

261  CONTINUE 

SOLVE  EQ  14 (SOLUTIONS) 

DO  271  J=1 , MM 
DO  269  L=l,  2 

269  XN( J, L)=XO( J, L)+DX( J, L) 

271  CONTINUE 

SOLVE  EQ  16 

(CALCULATE  THE  HIGHEST  DERIVATIVES) 

DO  281  1*1,2 
281  DW(I)*0.0 

DO  291  J*1,MM 
DO  285  K*l,2 
DV(K)*XN( J, K) 

DO  285  L*l,2 

285  DH(K, L)=A( J, K, L) 

CALL  MULV(DM,DV,DU, 2) 

DO  286  K*l,2 

286  DW(K)*DU(K)+DW(K) 

291  CONTINUE 

DO  301  K*l,2 

301  DW(K)=F(K, INDEX)+FDD(K, INDEX) -DW(K) 
CALL  MULV(AMI,DW,DV,2) 

DO  311  K*l,2 
311  XN(M,K)*DV(K) 


WRITE(6, * )XN(M, 1) ,XN(M, 2 ) 

PRINT  OUT  THE  DISPLACEMENT 

XRE*XN(1,2)-XN(1,1) 

WRITE ( 16,26 )T,XN( 1, 1 ) , XN( 1, 2 ) 
VEL1*XN(2, 1) 

VEL2*XN(2,2)-XN(2, 1) 
REACC*XN(3,2)-XN(3, 1) 

STORE  VALUES  AND  SET  XO  EQUAL  TO  XN 

ENED1 ( INDEX) =XN(2 , 1 ) * ( F( 2 , INDEX) - 
♦2 . 0*AM1*XN( 3 , 1 ) -AM2 *XN( 3 , 2 ) ) 

ENED2 ( INDEX)*XN( 2 , 2 ) * ( F( 2 , INDEX ) - 


+AM1*XN(3,1)-AM2*XN(3,2)) 

DO  321  J=1 , M 
DO  321  K=l,2 
321  X0( J,K)=XN( J,K) 

999  CONTINUE 

CALL  SIMP(ENED1,ENER1,DT,NT) 

CALL  S IMP ( ENED2 , ENER2 , DT , NT ) 
ENER=ENER1 +ENER2 
WRITE ( 6, 26)  ENER1 , ENER2 , ENER 
22  FORMAT (2E12. 4) 

26  FORMAT (3E12. 4) 

STOP 

END 

COMPUTE  THE  AREA  UNDER  THE  CURVE 

SUBROUTINE  SIMP(X, E, D, N) 
DIMENSION  X(1024) 

N2=N/2-l 

E=X(1) 

DO  1  I=1,N2 

IA=2*I 

IB=IA+1 


E=E+4.0*X(IA) 

E=E+2 . 0*X( IB) 

1  CONTINUE 
N3SN* 1 

E=E+0 . 5*X(N3 ) +1 . 5*X(N) 

E=D*E/3 . 0 
RETURN 
END 
C 

C  MATRIX  MULTIPLICATION 
C 

SUBROUTINE  MULT(A,B,C,N) 
DIMENSION  A(N,N),B(N,N),C(N,N) 
DO  1  1=1, N 
DO  1  J=1 , N 
C(I, J)*0.0 
DO  1  K=1,N 

C(I,J)=C(I,J)+A(I,K)*B(K,J) 

1  CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE  MULV(A,B,C,N) 
DIMENSION  A(N,N),B(N),C(N) 

DO  1  1=1, N 
C(I)=0.0 
DO  1  J=1,N 

C(I)=C(I)+A(I,J)*B(J) 


1  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  MATINV(C,N,D) 

C  MATRIX  INVERSION  C- INPUT  D- OUTPUT 
DIMENSION  C(N,N) ,D(N, N) 

DO  10  J=1,N 
DO  10  K=1,N 

10  D(J,K)sO.O 
DO  11  K=1 , N 

11  D(K,K)*1.0 
Pl*1.0 

DO  55  1=1, N 
P2=C(I, I) 

DO  40  J=1,N 
C(I, J)=C(I, J)/P2 
40  D(I,J)=D(I, J)/P2 
DO  51  IC=1 , N 
P3=-C( IC, I ) 

DO  50  K=1 , N 
IF( 10-1)21,51,21 
21  C(IC,K)=C(I,K)*P3+C(IC,K) 

50  D( IC, K)=D( I ,  K) *P3+D( IC, K) 

51  CONTINUE 
P1=P2*P1 

IF( ( I+2)-N)55, 53, 55 
53  DET=P1*( (C(I+1, l+l)*C(I+2, 1+2) )- 
+(C(I+2,I+1)*C(I+1, 1+2))) 

55  CONTINUE 

DO  70  IT=1,N 
DO  70  IS=1,N 
70  C(IT,IS)=D(IT,IS) 

RETURN 

END 


i 


184 


1 

8 


i 

i 

1 


k 


Q  ****************************************** 

C  *  SELECTED  SUBROUTINES  * 

C  ****************************************** 

C 

SUBROUTINE  SPCOEF(N,XN, FN, S, INDEX) 

DIMENSION  XN(N),FN(N),S(N),INDEX(N) 

DIMENSION  RHO( 1024) , TAU( 1024) 

NM1=N-1 

C 

C  ARRANGE  THE  NODES  XN  IN  INCREASING  ORDER. 

C  STORE  THE  ORDER  IN  THE  ARRAY  INDEX. 

C 

DO  1  1=1 , N 

1  INDEX ( I )=I 
DO  3  1=1 , NM1 
IP1=I+1 

DO  2  J=IP1,N 
11= INDEX ( I ) 

I J=INDEX( J) 

IF(XN( II ) . LE.XN( I J) )GO  TO  2 
ITEMP=INDEX( I ) 

INDEX ( I )=INDEX( J) 

INDEX( J)=ITEMP 

2  CONTINUE 

3  CONTINUE 
NM2=N-2 

C 

C  CALCULATE  THE  ELEMENTS  OF  THE  ARRAYS  RHO  AND  TAU. 
C 

RHO ( 2 ) =0 . 0 
TAU(2)=0.0 
DO  4  1=2 , NMl 
IIM1=INDEX( 1-1) 

I I = INDEX ( I ) 

IIP1=INDEX( 1+1) 

HIM1=XN( I I ) -XN( I IM1 ) 

HI=XN( IIP1)-XN( II ) 

TEMP= ( HIM1/HI ) * ( RHO ( I ) +2 . 0 ) +2 . 0 
RHO(I+l)=- 1.0/TEMP 
D=6.0*((FN(IIP1)-FN(II))/HI- 
+(FN( II )-FN( IIM1 ) )/HIMl )/HI 

4  TAU( I+1)=(D-HIM1*TAU( I )/HI )/TEMP 
C 

C  COMPUTE  ARRAY  OF  SECOND  DERIVATIVES  S  FOR 
C  THE  NATURAL  SPLINE 
C 

S(1)=0.0 
S{N)=0.0 
DO  5  1=1, NM2 
IB=N-I 

5  S( IB)*RH0( IB+1 ) *S( IB+1 ) +TAU( IB+1 ) 

RETURN 


R 


A  A 
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END 

FUNCTION  SPLINE(N,XN,FN,S, INDEX. X) 

DIMENSION  XN(N) , FN(N) , S(N) , INDEX(N) 

C 

C  IF  X.LT.XN( INDEX( 1) ) ,  APPROXIMATE  FUNCTION  BY 
C  THE  STRAIGHT  LINE  WHICH  PASSES  THROUGH  THE  POINT 
C  (XN( INDEX ( 1 ) ) , FN( IND£X( 1 ) ) )AND  WHOSE 

C  SLOPE  AGREES  WITH  THE  SPLINE  AT  THAT  POINT 
C 

11= INDEX ( 1 ) 

IF(X.GE.XN( II) )GO  TO  1 
I2=INDEX( 2 ) 

H1=XN( 12 )-XN( II) 

SPLINE=FN( I1)+(X-XN( 11))*((FN(I2)-FN(I1)) 
+/Hl-Hl*S(2)/6.0) 

RETURN 

IF  X.GE.XN( INDEX(N) ) ,  APPROXIMATE  FUNCTION  BY 
THE  STRAIGHT  LINE  WHICH  PASSES  THROUGH  THE  POINT 
(XN( INDEX(N) ) , FN( INDEX(N) ) )  AND  WHOSE  SLOPE  AGREES 
WITH  THE  SLOP  OF  THE  SPLINE  AT  THAT  POINT. 

1  IN=INDEX(N) 

IF(X. LE.XN( IN) )GO  TO  2 
INM1= INDEX (N-l) 

HNM1=XN( IN) -XN( INM1 ) 

SPLINE=FN ( IN ) + ( X-XN ( IN) )*( (FN( IN)-FN( INM1) ) 
+/HNM1+HNM1*S (N-l )/6 . 0) 

RETURN 

FOR  XN( INDEX ( 1 ) ) . LE . X . LE . XN ( INDEX ( N ) )  CALCULATE 
SPLINE  FIT. 

2  DO  3  1=2, N 
II=INDEX( I) 

IF(X. LE.XN( I I ) )G0  TO  4 

3  CONTINUE 

4  L=I-1 
IL*INDEX(L) 

ILP1»INDEX(L+1) 

A=XN( ILP1 ) -X 
B=X-XN{ IL) 

HL=XN( ILP1 ) -XN( IL) 

SPLINE=A*S( L) * ( A**2/HL-HL)/6 . 0+B*S ( L+l ) * 
+(B**2/HL-HL)/6.0-(A*FN( IL)+B*FN( ILP1) )/HL 
RETURN 
END 

FFT  PROGRAM 

******************************************** 
MULTIPLIED  BY  T  WHEN  USING  THE  FORWARD  FFT 
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C  DIVIDED  BY  DT  FOR  BACKWARD  FFT 
C  SIGN  =-l.  FOR  DFT,  SIGN=1.  FOR  ID FT 

C  ******************************************** 

C 

SUBROUTINE  FFT1(A,N,NB, SIGN) 

COMPLEX  A(NB),U,W,T 

DIVIDE  ALL  ELEMENT  BY  NB 

DO  1  J=1 , NB 

1  A(J)=A(J)/NB 

REORDER  SEQUENCE 

NBD2=NB/2 
NBM1=NB-1 
J=1 

DO  4  L=1 , NBM1 
IF(L.GE.J)  GO  TO  2 
T=A( J) 

A(J)=A(L) 

A(L)=T 

2  K=NBD2 

3  IF(K.GE.J)  GO  TO  4 
J=J-K 
'K—K/2 
GO  TO  3 

4  J=J+K 
CALCULATE  FFT 

PI=3 . 141592653589793 
DO  6  M=l,  N 
U=(l. 0,0.0) 

ME=2**M 
K=ME/2 

W=CMPLX( COS ( PI/K ) , S IGN*SIN( P I/K ) ) 

DO  6  J=1 , K 
DO  5  L= J ,  NB ,  ME 
LPK=L+K 
T*A(LPK)*U 
A(LPK)=A(L)-T 

5  A(L)—A(L) +T 

6  U=U*W 
RETURN 
END 

SUBROUTINE  SIMP (QF, A, B, ACC, ANS, ERROR, AREA, I FLAG) 

C 

C  SIMP  IS  AN  ADAPTIVE,  ITERATIVE  CODE  BASED 
C  ON  SIMPSON'S  RULE.  IT  IS  DESIGNED  TO  EVALUATE  THE 
C  DEFINITE  INTEGRAL  OF  A  CONTINUOUS  FUNCTION  WITH 
C  FINITE  LIMITS  OF  INTEGRATION. 


_ _ •  ^  JL—r  _  m  . _ * _ -  »  ■  m  ■ 


F  -  NAME  OF  FUNCTION  WHOSE  INTEGRAL  IS  DESIRED. 

THE  FUNCTION  NAME  F  MUST  APPEAR  IN  AN  EXTERNAL 
STATEMENT  IN  THE  CALLING  PROGRAM. 

A,B  -  LOWER  AND  UPPER  LIMITS  OF  INTEGRATION. 

ANS-  APPROXIMATE  VALUE  OF  THE  INTEGRAL  OF  F(X) 

FROM  A  TO  B. 

AREA  -  APPROXIMATE  VALUE  OF  THE  INTEGRAL  OF 
ABS(F(X) )  FROM  A  TO  B. 

ERROR  -  ESTIMATED  ERROR  OF  ANS.  USER  MAY  WISH 
TO  EXTRAPOLATE  BY  FORMING  ANS+ERROR  TO  GET  WHAT  IS 
OFTEN  A  MORE  ACCURATE  RESULT,  BUT  NOT  ALWAYS. 

ACC  -  DESIRED  ACCURACY  OF  ANS.  CODE  TRIES  TO  MAKE 
ABS ( ERROR ) . LE . ACC*ABS ( AREA ) . 

I FLAG  =  1  FOR  NORMAL  RETURN. 

=  2  IF  IT  IS  NECESSARY  TO  GO  TO  30  LEVELS  OR 
USE  A  SUB INTERVAL  TOO  SMALL  FOR  MACHINE  WORD 
LENGTH.  ERROR  MAY  BE  UNRELIABLE  IN  THIS  CASE 
=  3  IF  MORE  THAT  2000  FUNCTION  EVALUATIONS 
ARE  USED.  ROUGH  APPROXIMATIONS  ARE  USED  TO 
COMPLETE  THE  COMPUTATIONS  AND  ERROR  IS 
USUALLY  UNRELIABLE. 

DIMENSION  FV(5) ,L0RR(30) , FIT (30) ,F2T(30) ,F3T(30) 
DIMENSION  ARESTT( 30) , ESTT ( 30 ) , EPST ( 30 ) , PSUM( 30 ) 
DIMENSION  DAT (30) 

SET  U  TO  APPROXIMATELY  THE  UNIT  ROUND-OFF  OF 
SPECIFIC  MACHINE  (HERE  IBM  360/67) 

U  =  9.0E-7 

INITIALIZE 

FOURU=4.0*U 
IFLAG=1 
EPS=ACC 
ERROR=0 . 0 
LVL=1 

L0RR(LVL)=1 
PSUM(LVL)=0 . 0 
ALPHA=A 
DA=B-A 
AREA—0.0 
AREST=0 . 0 
FV(1)=QF( ALPHA) 

FV ( 3 ) =QF ( ALPHA+O . 5*DA) 

FV ( 5 ) *QF ( ALPHA+DA ) 

K0UNT=3 

wt=d;  /6 . 0 

EST**  :*(FV(1)+4.0*FV(3)+FV(5)) 
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C  'BASIC  STEP'.  HAVE  ESTIMATE  EST  OF  INTEGRAL 
C  ON  ( ALPHA , ALPHA+DA ) .  BISECT  AND  COMPUTE 
C  ESTIMATES  ON  LEFT  AND  RIGHT  HALF  INTERVALS. 

C  SIMILARLY  TREAT  INTEGRAL  OF  ABS(F(I>).  SUM  IS 
C  BETTER  VALUE  FOR  INTEGRAL  AND  DIFF/15.0  IS 
C  APPROXIMATELY  ITS  ERROR. 

C 

1  DX=0 . 5  *DA 

FV ( 2 ) =QF ( ALPHA+O . 5 *DX ) 

FV ( 4 ) *QF ( ALPHA+ 1 . 5 *DX ) 

KOUNT=KOUNT+2 
WT— DX/6 . 0 

E STL=WT * ( FV ( 1 ) + 4 . 0 * FV ( 2 ) + FV ( 3 ) ) 

ESTR=WT* ( FV( 3 ) +4 . 0*FV( 4 ) +FV ( 5 ) ) 

SUM=ESTL+ESTR 

ARESTL=WT*(ABS(FV( 1 ) )+ABS (4. 0*FV(2 ) )+ABS(FV(3 ) ) ) 
ARESTR=WT* ( ABS ( FV( 3 ) ) +ABS ( 4 . 0*FV ( 4 ) )+ABS(FV(5) ) ) 
AREA=AREA+ ( ( ARESTL+ARESTR ) - AREST ) 

DIFF=EST-SUM 

IF  ERROR  IS  ACCEPTABLE,  GO  TO  2.  IF  INTERVAL 
IS  TOO  SMALL  OR  TOO  MANY  LEVELS  OR  TOO  MANY 
FUNCTION  EVALUATIONS,  SET  A  FLAG 
AND  GO  TO  2  ANYWAY. 

IF(ABS(DIFF) .LE. EPS* ABS (AREA) )GO  TO  2 
IF(ABS(DX) .LE.FOURU* ABS (ALPHA) )GO  TO  5 
IF(LVL.GE.30)GO  TO  5 
IF(KOUNT.GE.2000)GO  TO  6 

HERE  TO  RAISE  LEVEL.  STORE  INFORMATION  TO 
PROCESS  RIGHT  HALF  INTERVAL  LATER.  INITIALIZE  FOR 
'BASIC  STEP'  SO  AS  TO  TREAT  LEFT  HALF  INTERVAL. 

LVL=LVL+1 
LORR ( LVL ) =0 
F1T(LVL)=FV(3 ) 

F2T ( LVL ) =FV ( 4 ) 

F3T( LVL)=FV( 5 ) 

DA=DX 

DAT(LVL)=DX 
AREST=ARESTL 
ARESTT ( LVL ) =ARESTR 
EST=ESTL 
ESTT ( LVL ) =ESTR 
EPS*EPS/1 . 4 
EPST ( LVL ) =EPS 
FV(5 )=FV(3 ) 

FV ( 3 )=FV ( 2 ) 

GO  TO  1 

ACCEPT  APPROXIMATE  INTEGRAL  SUM.  IF  IT  WAS  ON 


A  LEFT  INTERVAL  GO  TO  'MOVE  RIGHT' .  IF  A  RIGHT 
INTERVAL,  ADD  RESULTS  TO  FINISH  AT  THIS  LEVEL. 
ARRAY  LORR  (MNEMONI"  FOR  LEFT  OR  RIGHT  TELLS 
WHETHER  LEFT  OR  RIGHT  INTERVAL  AT  EACH  LEVEL. 

ERROR=ERROR+D IFF/15 . 0 
IF ( LORR ( LVL ) . EQ . 0 ) GO  TO  4 
SUM=PSUM ( LVL ) +  SUM 
LVL=LVL-1 

IF ( LVL . GT . 1 ) GO  TO  3 

ANS=SUM 

RETURN 

'MOVE  RIGHT'.  RESTORE  SAVED  INFORMATION  TO 
PROCESS  RIGHT  HALF  INTERVAL. 

PSUM ( LVL ) =SUM 
LORR ( LVL )=1 
ALPHA= ALPHA + DA 
DA=DAT ( LVL ) 

FV ( 1 ) =F1T ( LVL ) 

FV ( 3 ) =F2T ( LVL ) 

FV(5)=F3T(LVL) 

AREST=ARESTT ( LVL ) 

EST=ESTT ( LVL ) 

EPS=EPST ( LVL ) 

GO  TO  1 

ACCEPT  'POOR'  VALUE.  SET  APPROPRIATE  FLAGS. 

IFLAG=2 
GO  TO  2 
IFLAG=3 
GO  TO  2 
END 


SUBROUTINE  SPCOEF(N,XN, FN, S, INDEX) 
DIMENSION  XN (N),FN(N),S(N), INDEX ( N ) 
DIMENSION  RHO( 1024) ,TAU( 1024) 

NM1=N- 1 

ARRANGE  THE  NODES  XN  IN  INCREASING  ORDER. 
STORE  THE  ORDER  IN  THE  ARRAY  INDEX. 

DO  1  1=1, N 
INDEX ( I )=I 
DO  3  1=1, NM1 
IP1=I+1 
DO  2  J=IP1,N 
II=INDEX( I ) 

I J=INDEX( J) 
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IF(XN( II ) . LE.XN( IJ> )GO  TO  2 
ITEMP=INDEX( I) 

INDEX( I )=INDEX( J) 

INDEX( J)=ITEMP 

2  CONTINUE 

3  CONTINUE 
NM2-N-2 

CALCULATE  THE  ELEMENTS  OF  THE  ARRAYS 
RHO  AND  TAU. 

RHO ( 2 ) =0 . 0 
TAU ( 2 ) =0 . 0 
DO  4  1=2, NM1 
I IM1=INDEX( I -1 ) 

I I = INDEX ( I ) 

IIP1=INDEX( 1+1) 

HIM1=XN( I I ) -XN( I IM1 ) 

HI=XN( IIP1)-XN( II ) 

TEMP=(HIM1/HI ) * ( RHO ( I ) +2 . 0 ) +2 . 0 

RHO ( I +1)=-1. 0/TEMP 
D=6.0*( (FN( IIP1)-FN( II ) )/HI 
+-(FN( II)-FN{ IIM1) )/HIMl)/HI 

4  TAU(I+1)=(D-HIM1*TAU( I)/HI)/TEMP 

COMPUTE  ARRAY  OF  SECOND  DERIVATIVES  S  FOR 
THE  NATURAL  SPLINE. 

S(1)=0.0 
S ( N ) =0 . 0 
DO  5  1=1 , NM2 
IB=N- I 

5  S(IB)=RHO(IB+l)*S(IB+l)+TAU(IB+l) 

RETURN 
END 

FUNCTION  SPLINE(N,XN,FN,S, INDEX, X) 

DIMENSION  XN(N) , FN(N) , S (N) , INDEX(N) 

C 

C  IF  X. LT.XN( INDEX ( 1 ) ) ,  APPROXIMATE  FUNCTION 

C  BY  THE  STRAIGHT  LINE  WHICH  PASSES  THROUGH  THE 
C  POINT  ( XN ( INDEX ( 1 ) ) , FN ( INDEX ( 1 ) ) )  AND  WHOSE  SLOPE 
C  AGREES  WITH  THE  SLOPE  OF  THE  SPLINE  AT  THAT  POINT. 
C 

I 1=INDEX( 1 ) 

IF(X.GE.XN( II) )GO  TO  1 
I2=INDEX( 2 ) 

H1=XN( 12 )-XN( II) 

SPLINE=FN(I1)+(X-XN(I1) )*( (FN( 12 ) -FN( II ) ) 
+/H1-H1*S(2 )/6. 0) 

RETURN 

C  IF  X . GE . XN ( INDEX ( N ) ) ,  APPROXIMATE  FUNCTION  BY 
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THE  STRAIGHT  LINE  WHICH  PASSES  THROUGH  THE  POINT 
(XN( INDEX(N) ) , FN( INDEX(N) ) )  AND  WHOSE  SLOPE 
AGREES  WITH  THE  SLOPE  OF  THE  SPLINE  AT  THAT  POINT. 

1  IN=INDEX(N) 

IF(X. LE . XN( IN) )GO  TO  2 
INM1= INDEX (N-l) 

HNM1=XN( IN) -XN( INM1 ) 

SPLINE=FN(IN)+(X-XN(IN) )*( (FN(IN)-FN(INMl) )/HNMl 
++HNMl*S(N-l)/6 . 0 ) 

RETURN 

FOR  XN( INDEX ( 1 ) ) . LE .X. LE .XN( INDEX (N) ) 

CALCULATE  SPLINE  FIT. 

2  DO  3  1=2 , N 
I I=INDEX( I ) 

IF(X.LE.XN( II ) )GO  TO  4 

3  CONTINUE 

4  L=I-1 
IL=INDEX(L) 

ILP1=INDEX( L+l ) 

A=XN( ILP1 ) -X 
B=X-XN( IL) 

HL=XN( ILP1)-XN( IL) 

SPLINE=A*S(L) * (A**2/HL-HL)/6. 0+B*S(L+l) * (B**2/HL-HL) 
+/6.0-(A*FN( IL)+B*FN( ILP1) )/HL 
RETURN 
END 

THIS  SUBROUTINE  IS  FOR  N*N  MATRIX  INVERSION 

SUBROUTINE  INVERS(B,N, A) 

DOUBLE  PRECISION  A(2,2) ,B(2,2) ,C(2,2) 

DOUBLE  PRECISION  AMAX, TEMP, PIVOT 
DIMENSION  INDEX(4, 2 ) 

IF(N.GT. 40)  GO  TO  134 
DO  90  I®1, N 
DO  90  J=1,N 
A(I,J)=B(I,J) 

90  CONTINUE 

DO  107  1=1, N 
DO  107  J=1 ,  N 

107  B( I, J)=A( I , J) 

DO  108  1=1, N 

108  INDEX ( I , 1 )=0 
11=0 

109  AMAX=-1 . 0D0 
DO  110  1=1, N 

IF( INDEX ( 1,1))  110,111,110 
111  DO  112  J=1,N 
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I F ( INDEX( J, 1) )  112,113,112 

113  TEMP=DABS(A(I, J) ) 
IF(TEMP-AMAX)  112,112,114 

114  IROW=I 
ICOL=J 
AMAX=TEMP 

112  CONTINUE 
110  CONTINUE 

IF(AMAX)  225,115,116 
116  INDEX ( ICOL,  1)=IR0W. 

IF( IROW-ICOL)  119,118,119 

119  DO  120  J=1,N 
TEMP=A(  IROW,  J) 

A( IROW, J)=A( ICOL, J) 

120  A( ICOL, J)=TEMP 
11=11+1 

INDEX( II ,2 )=ICOL 
118  PIVOT=A( ICOL, ICOL) 

A( ICOL, IC0L)=1 . ODO 
PIV0T=1 . ODO/PIVOT 
DO  121  J=1 , N 

121  A( ICOL, J)=A( ICOL, J)*PIVOT 
DO  122  1=1, N 
IF(I-ICOL)  123,122,123 

123  TEMP=A(I, ICOL) 

A( I , ICOL ) =0 . ODO 
DO  124  J=1,N 

124  A( I , J)=A( I , J)-A( ICOL, J) *TEMP 

122  CONTINUE 
GO  TO  109 

125  IC0L=INDEX(II,2) 

IROW=INDEX( ICOL, 1 ) 

DO  126  1=1, N 
TEMP=A( I , IROW) 

A( I, IROW)=A( I , ICOL) 

126  A( I , ICOL)=TEMP 
11=11-1 

225  IF(II)  125,127,125 

127  DO  130  1=1, N 
DO  130  J=1,N 
C(I, J)=0.0D0 
DO  130  K=1 , N 

130  C(I,J)=C(I, J)+B(I,K)*A(K, J) 
GO  TO  134 

115  WRITE( 6, 133 ) 

133  FORMAT ('  ZERO  PIVOT') 

134  RETURN 
END 


LEAST  SQUARE  POLYNOMINAL  FITTING 


non 
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C  Y  :  INPUT 

C  MP  :  DEGREE  OF  POLYNOMINAL 

C  N  : NUMBER  OF  POINT  FOR  THE  FITTING 

C  XI  :  THE  LOW  X-AXIS  VALUE 

C  X2  :  THE  UPPER  X-AXIS  VALUE 

C  A( I ) :  THE  COEFFICIENT  OF  THE  POLYNOMINAL 

C  A( 1) :  THE  CONSTANT  COEFFICIENT 

C 

SUBROUTINE  FIT1(Y,M?,N,X1,X2, A) 

DIMENSION  XX(20,3),Y(20),X(20),XP(3,3),A(3) 
DIMENSION  EPS (20),D(3,3),B(3) 

DX= (X2-X1 ) /FLOAT ( N- 1 ) 

DO  10  1=1 , N 

10  X(I)=X1+FL0AT(I-1)*DX 
DO  1  1*1, N 
XX(I,1)*1.0 
DO  1  J=2,MP 

1  XX(I, J)=X(I)**(J-1) 

DO  2  1=1 ,MP 

DO  2  J*I,MP 
XP(I, J)=0.0 
DO  2  K=1 , N 

2  XP(I,J)=XP(I, J)+XX(K,I)*XX(K, J) 

DO  3  1=2, MP 

IM=I-1 
DO  3  J=l, IM 

3  XP(I,J)*XP(J, I) 

CALL  MATINV(XP,MP,D) 

DO  4  1=1, MP 

B( I )=0. 0 
DO  4  J=1,N 

4  B( I )=B( I )+XX( J, I ) *Y( J) 

DO  5  1=1, MP 

A( I )=0. 0 
DO  5  J=1 , MP 

5  A(I)=A(I)+D(I, J)*B(J) 

RETURN 
END 

MATRIX  INVERSION 

SUBROUTINE  MATINV(C,N,D) 

C  MATRIX  INVERSION  C- INPUT  D- OUTPUT 
DIMENSION  C(N,N) ,D(N,N) 

DO  10  J=1,N 
DO  10  K=1,N 

10  D(J,K)=0.0 
DO  11  K=1,N 

11  D(K, K)=l . 0 
Pl=1.0 

DO  55  1=1, N 
P2=C( 1,1) 


DO  40  J=1 ,  N 
C(I, J)=C(I, J)/P2 
40  D(I,J)=D(I,J)/P2 
DO  51  IC=1,N 
P3=-C( IC, X ) 

DO  50  K=1 , N 
IF( 10-1)21/51,21 
21  C(IC,K)=C(I,K)*P3+C(IC,K) 

50  D( IC, K)=D( I ,K) *P3+D( IC, K) 

51  CONTINUE 
P1*P2*P1 

IF((I+2)-N)55,53,55 
53  DET=Pl*((C(I+l,I+l)*C(I+2,I+2)) 
+(C(I+2,I+l)*C(I+l,I+2))) 

55  CONTINUE 

DO  70  IT=1,N 
DO  70  IS=1,N 
70  C( IT, IS)*D( IT, IS) 

RETURN 


